Context Navigation

readsulfate.F @ 923

Last change on this file since 923 was 804, checked in by Laurent Fairhead, 17 years ago
Un probleme sur le calcul des niveaux verticaux L.Li LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.8 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
[766]	4	SUBROUTINE readsulfate (r_day, first, sulfate_p)
[776]	5	USE dimphy, ONLY : klev
	6	USE mod_grid_phy_lmdz, klon=>klon_glo
	7	USE mod_phys_lmdz_para
[524]	8	IMPLICIT none
	9
	10	c Content:
	11	c --------
	12	c This routine reads in monthly mean values of sulfate aerosols and
	13	c returns a linearly interpolated dayly-mean field.
	14	c
	15	c
	16	c Author:
	17	c -------
	18	c Johannes Quaas (quaas@lmd.jussieu.fr)
	19	c 26/04/01
	20	c
	21	c Modifications:
	22	c --------------
	23	c 21/06/01: Make integrations of more than one year possible ;-)
	24	c ATTENTION!! runs are supposed to start with Jan, 1. 1930
	25	c (rday=1)
	26	c
	27	c 27/06/01: Correction: The model always has 360 days per year!
	28	c 27/06/01: SO4 concentration rather than mixing ratio
	29	c 27/06/01: 10yr-mean-values to interpolate
	30	c 20/08/01: Correct the error through integer-values in interpolations
	31	c 21/08/01: Introduce flag to read in just one decade
	32	c
	33	c Include-files:
	34	c --------------
	35	#include "YOMCST.h"
	36	#include "chem.h"
	37	#include "dimensions.h"
[766]	38	cym#include "dimphy.h"
[524]	39	#include "temps.h"
	40	c
	41	c Input:
	42	c ------
[782]	43	REAL r_day ! Day of integration
[524]	44	LOGICAL first ! First timestep
	45	! (and therefore initialization necessary)
	46	c
	47	c Output:
	48	c -------
[782]	49	REAL sulfate_p(klon_omp,klev)
	50	REAL sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
[524]	51	! from file) [ug SO4/m3]
	52	c
	53	c Local Variables:
	54	c ----------------
	55	INTEGER i, ig, k, it
[640]	56	INTEGER j, iday, ny, iyr, iyr1, iyr2
[524]	57	parameter (ny=jjm+1)
	58
[640]	59	CJLD INTEGER idec1, idec2 ! The two decadal data read ini
[524]	60	CHARACTER*4 cyear
	61
	62	INTEGER im, day1, day2, im2
[782]	63	REAL so4_1(iim, jjm+1, klev, 12)
	64	REAL so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions
[524]	65
[782]	66	REAL, allocatable,save :: so4(:, :, :) ! SO4 in right dimension
	67	REAL, allocatable,save :: so4_out(:, :)
[766]	68	c$OMP THREADPRIVATE(so4,so4_out)
[524]	69
	70	LOGICAL lnewday
	71	LOGICAL lonlyone
	72	PARAMETER (lonlyone=.FALSE.)
[766]	73	logical,save :: first2=.true.
	74	c$OMP THREADPRIVATE(first2)
[524]	75
[766]	76	c$OMP MASTER
	77	if (first2) then
	78
	79	allocate( so4(klon, klev, 12) )
	80	allocate( so4_out(klon, klev))
	81	first2=.false.
	82
	83	endif
	84
[776]	85	if (is_mpi_root) then
[766]	86
[524]	87	iday = INT(r_day)
	88
	89	! Get the year of the run
	90	iyr = iday/360
	91
	92	! Get the day of the actual year:
	93	iday = iday-iyr*360
	94
	95	! 0.02 is about 0.5/24, namly less than half an hour
	96	lnewday = (r_day-FLOAT(iday).LT.0.02)
	97
	98	! ---------------------------------------------
	99	! All has to be done only, if a new day begins!
	100	! ---------------------------------------------
	101
	102	IF (lnewday.OR.first) THEN
	103
	104	im = iday/30 +1 ! the actual month
	105	! annee_ref is the initial year (defined in temps.h)
	106	iyr = iyr + annee_ref
	107
	108	! Do I have to read new data? (Is this the first day of a year?)
	109	IF (first.OR.iday.EQ.1.) THEN
	110	! Initialize values
	111	DO it=1,12
	112	DO k=1,klev
	113	DO i=1,klon
	114	so4(i,k,it)=0.
	115	ENDDO
	116	ENDDO
	117	ENDDO
	118
[640]	119
	120	IF (iyr .lt. 1850) THEN
	121	cyear='.nat'
	122	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
	123	CALL getso4fromfile(cyear, so4_1)
	124	ELSE IF (iyr .ge. 2100) THEN
	125	cyear='2100'
	126	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
	127	CALL getso4fromfile(cyear, so4_1)
	128	ELSE
	129
	130	! Read in data:
[524]	131	! a) from actual 10-yr-period
	132
[640]	133	IF (iyr.LT.1900) THEN
	134	iyr1 = 1850
	135	iyr2 = 1900
	136	ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
	137	iyr1 = 1900
	138	iyr2 = 1920
	139	ELSE
	140	iyr1 = INT(iyr/10)*10
	141	iyr2 = INT(1+iyr/10)*10
[524]	142	ENDIF
[640]	143	WRITE(cyear,'(I4)') iyr1
	144	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
[524]	145	CALL getso4fromfile(cyear, so4_1)
	146
	147
	148	! If to read two decades:
	149	IF (.NOT.lonlyone) THEN
	150
	151	! b) from the next following one
[640]	152	WRITE(cyear,'(I4)') iyr2
	153	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
	154	CALL getso4fromfile(cyear, so4_2)
	155
[524]	156	ENDIF
[640]	157
[524]	158	! Interpolate linarily to the actual year:
	159	DO it=1,12
	160	DO k=1,klev
	161	DO j=1,jjm
	162	DO i=1,iim
	163	so4_1(i,j,k,it)=so4_1(i,j,k,it)
[640]	164	. - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
[524]	165	. * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
	166	ENDDO
	167	ENDDO
	168	ENDDO
	169	ENDDO
	170
	171	ENDIF !lonlyone
	172
	173	! Transform the horizontal 2D-field into the physics-field
	174	! (Also the levels and the latitudes have to be inversed)
	175
	176	DO it=1,12
	177	DO k=1, klev
	178	! a) at the poles, use the zonal mean:
	179	DO i=1,iim
	180	! North pole
	181	so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
	182	! South pole
	183	so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
	184	ENDDO
	185	so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
	186	so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
	187
	188	! b) the values between the poles:
	189	ig=1
	190	DO j=2,jjm
	191	DO i=1,iim
	192	ig=ig+1
	193	if (ig.gt.klon) write (,) 'shit'
[804]	194	so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	195	ENDDO
	196	ENDDO
	197	IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
	198	ENDDO ! Loop over k (vertical)
	199	ENDDO ! Loop over it (months)
	200
	201
	202	ENDIF ! Had to read new data?
	203
	204
	205	! Interpolate to actual day:
	206	IF (iday.LT.im*30-15) THEN
	207	! in the first half of the month use month before and actual month
	208	im2=im-1
	209	day2 = im2*30-15
	210	day1 = im2*30+15
	211	IF (im2.LE.0) THEN
	212	! the month is january, thus the month before december
	213	im2=12
	214	ENDIF
	215	DO k=1,klev
	216	DO i=1,klon
	217	sulfate(i,k) = so4(i,k,im2)
	218	. - FLOAT(iday-day2)/FLOAT(day1-day2)
	219	. * (so4(i,k,im2) - so4(i,k,im))
	220	IF (sulfate(i,k).LT.0.) THEN
	221	IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
	222	IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
	223	. write(,) 'so4(i,k,im2) - so4(i,k,im)',
	224	. so4(i,k,im2) - so4(i,k,im)
	225	IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
	226	stop 'sulfate'
	227	endif
	228	ENDDO
	229	ENDDO
	230	ELSE
	231	! the second half of the month
	232	im2=im+1
	233	IF (im2.GT.12) THEN
	234	! the month is december, the following thus january
	235	im2=1
	236	ENDIF
	237	day2 = im*30-15
	238	day1 = im*30+15
	239	DO k=1,klev
	240	DO i=1,klon
	241	sulfate(i,k) = so4(i,k,im2)
	242	. - FLOAT(iday-day2)/FLOAT(day1-day2)
	243	. * (so4(i,k,im2) - so4(i,k,im))
	244	IF (sulfate(i,k).LT.0.) THEN
	245	IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
	246	IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
	247	. write(,) 'so4(i,k,im2) - so4(i,k,im)',
	248	. so4(i,k,im2) - so4(i,k,im)
	249	IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
	250	stop 'sulfate'
	251	endif
	252	ENDDO
	253	ENDDO
	254	ENDIF
	255
	256
[640]	257	CJLD ! The sulfate concentration [molec cm-3] is read in.
	258	CJLD ! Convert it into mass [ug SO4/m3]
	259	CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
	260	! The sulfate mass [ug SO4/m3] is read in.
[524]	261	DO k=1,klev
	262	DO i=1,klon
[640]	263	CJLD sulfate(i,k) = sulfate(i,k)*masse_so4
	264	CJLD . /n_avogadro*1.e+12
[524]	265	so4_out(i,k) = sulfate(i,k)
	266	IF (so4_out(i,k).LT.0)
	267	. stop 'WAS SOLL DER SCHEISS ? '
	268	ENDDO
	269	ENDDO
	270	ELSE ! if no new day, use old data:
	271	DO k=1,klev
	272	DO i=1,klon
	273	sulfate(i,k) = so4_out(i,k)
	274	IF (so4_out(i,k).LT.0)
	275	. stop 'WAS SOLL DER SCHEISS ? '
	276	ENDDO
	277	ENDDO
	278
	279
	280	ENDIF ! Did I have to do anything (was it a new day?)
[766]	281
	282	endif ! phy_rank==0
[524]	283
[766]	284	c$OMP END MASTER
[782]	285	call Scatter(sulfate,sulfate_p)
[766]	286
[524]	287	RETURN
	288	END
	289
	290
	291
	292
	293
	294	c-----------------------------------------------------------------------------
	295	c Read in /calculate pre-industrial values of sulfate
	296	c-----------------------------------------------------------------------------
	297
[766]	298	SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p)
[776]	299	USE dimphy, ONLY : klev
	300	USE mod_grid_phy_lmdz, klon=>klon_glo
	301	USE mod_phys_lmdz_para
[524]	302	IMPLICIT none
	303
	304	c Content:
	305	c --------
	306	c This routine reads in monthly mean values of sulfate aerosols and
	307	c returns a linearly interpolated dayly-mean field.
	308	c
	309	c It does so for the preindustriel values of the sulfate, to a large part
	310	c analogous to the routine readsulfate above.
	311	c
	312	c Only Pb: Variables must be saved and don t have to be overwritten!
	313	c
	314	c Author:
	315	c -------
	316	c Johannes Quaas (quaas@lmd.jussieu.fr)
	317	c 26/06/01
	318	c
	319	c Modifications:
	320	c --------------
	321	c see above
	322	c
	323	c Include-files:
	324	c --------------
	325	#include "YOMCST.h"
	326	#include "chem.h"
	327	#include "dimensions.h"
	328	#include "temps.h"
	329	c
	330	c Input:
	331	c ------
[782]	332	REAL r_day ! Day of integration
[524]	333	LOGICAL first ! First timestep
	334	! (and therefore initialization necessary)
	335	c
	336	c Output:
	337	c -------
[782]	338	REAL pi_sulfate_p (klon_omp, klev)
[766]	339
[782]	340	REAL pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
[766]	341	! from fil
[524]	342	c
	343	c Local Variables:
	344	c ----------------
	345	INTEGER i, ig, k, it
	346	INTEGER j, iday, ny, iyr
	347	parameter (ny=jjm+1)
	348
[782]	349	INTEGER im, day1, day2, im2
	350	REAL pi_so4_1(iim, jjm+1, klev, 12)
[766]	351
[782]	352	REAL, allocatable,save :: pi_so4(:, :, :) ! SO4 in right dimension
	353	REAL, allocatable,save :: pi_so4_out(:, :)
[766]	354	c$OMP THREADPRIVATE(pi_so4,pi_so4_out)
[524]	355
	356	CHARACTER*4 cyear
	357	LOGICAL lnewday
[766]	358	logical,save :: first2=.true.
	359	c$OMP THREADPRIVATE(first2)
[524]	360
[766]	361	c$OMP MASTER
	362	if (first2) then
	363
	364	allocate( pi_so4(klon, klev, 12) )
	365	allocate( pi_so4_out(klon, klev))
	366	first2=.false.
	367
	368	endif
	369
[776]	370	if (is_mpi_root) then
[766]	371
[524]	372
	373
	374	iday = INT(r_day)
	375
	376	! Get the year of the run
	377	iyr = iday/360
	378
	379	! Get the day of the actual year:
	380	iday = iday-iyr*360
	381
	382	! 0.02 is about 0.5/24, namly less than half an hour
	383	lnewday = (r_day-FLOAT(iday).LT.0.02)
	384
	385	! ---------------------------------------------
	386	! All has to be done only, if a new day begins!
	387	! ---------------------------------------------
	388
	389	IF (lnewday.OR.first) THEN
	390
	391	im = iday/30 +1 ! the actual month
	392
	393	! annee_ref is the initial year (defined in temps.h)
	394	iyr = iyr + annee_ref
	395
	396
	397	IF (first) THEN
	398	cyear='.nat'
	399	CALL getso4fromfile(cyear,pi_so4_1)
	400
	401	! Transform the horizontal 2D-field into the physics-field
	402	! (Also the levels and the latitudes have to be inversed)
	403
	404	! Initialize field
	405	DO it=1,12
	406	DO k=1,klev
	407	DO i=1,klon
	408	pi_so4(i,k,it)=0.
	409	ENDDO
	410	ENDDO
	411	ENDDO
	412
	413	write (,) 'preind: finished reading', FLOAT(iim)
	414	DO it=1,12
	415	DO k=1, klev
	416	! a) at the poles, use the zonal mean:
	417	DO i=1,iim
	418	! North pole
	419	pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
	420	! South pole
	421	pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
	422	ENDDO
	423	pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
	424	pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
	425
	426	! b) the values between the poles:
	427	ig=1
	428	DO j=2,jjm
	429	DO i=1,iim
	430	ig=ig+1
	431	if (ig.gt.klon) write (,) 'shit'
[804]	432	pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	433	ENDDO
	434	ENDDO
	435	IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
	436	ENDDO ! Loop over k (vertical)
	437	ENDDO ! Loop over it (months)
	438
	439	ENDIF ! Had to read new data?
	440
	441
	442	! Interpolate to actual day:
	443	IF (iday.LT.im*30-15) THEN
	444	! in the first half of the month use month before and actual month
	445	im2=im-1
	446	day1 = im2*30+15
	447	day2 = im2*30-15
	448	IF (im2.LE.0) THEN
	449	! the month is january, thus the month before december
	450	im2=12
	451	ENDIF
	452	DO k=1,klev
	453	DO i=1,klon
	454	pi_sulfate(i,k) = pi_so4(i,k,im2)
	455	. - FLOAT(iday-day2)/FLOAT(day1-day2)
	456	. * (pi_so4(i,k,im2) - pi_so4(i,k,im))
	457	IF (pi_sulfate(i,k).LT.0.) THEN
	458	IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
	459	IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
	460	. write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
	461	. pi_so4(i,k,im2) - pi_so4(i,k,im)
	462	IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
	463	stop 'pi_sulfate'
	464	endif
	465	ENDDO
	466	ENDDO
	467	ELSE
	468	! the second half of the month
	469	im2=im+1
	470	day1 = im*30+15
	471	IF (im2.GT.12) THEN
	472	! the month is december, the following thus january
	473	im2=1
	474	ENDIF
	475	day2 = im*30-15
	476
	477	DO k=1,klev
	478	DO i=1,klon
	479	pi_sulfate(i,k) = pi_so4(i,k,im2)
	480	. - FLOAT(iday-day2)/FLOAT(day1-day2)
	481	. * (pi_so4(i,k,im2) - pi_so4(i,k,im))
	482	IF (pi_sulfate(i,k).LT.0.) THEN
	483	IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
	484	IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
	485	. write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
	486	. pi_so4(i,k,im2) - pi_so4(i,k,im)
	487	IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
	488	stop 'pi_sulfate'
	489	endif
	490	ENDDO
	491	ENDDO
	492	ENDIF
	493
	494
[640]	495	CJLD ! The sulfate concentration [molec cm-3] is read in.
	496	CJLD ! Convert it into mass [ug SO4/m3]
	497	CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
[524]	498	DO k=1,klev
	499	DO i=1,klon
[640]	500	CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
	501	CJLD . /n_avogadro*1.e+12
[524]	502	pi_so4_out(i,k) = pi_sulfate(i,k)
	503	ENDDO
	504	ENDDO
	505
	506	ELSE ! If no new day, use old data:
	507	DO k=1,klev
	508	DO i=1,klon
	509	pi_sulfate(i,k) = pi_so4_out(i,k)
	510	ENDDO
	511	ENDDO
	512
	513
	514	ENDIF ! Was this the beginning of a new day?
[766]	515
[776]	516	endif ! is_mpi_root==0
[766]	517
	518	c$OMP END MASTER
[782]	519	call Scatter(pi_sulfate,pi_sulfate_p)
[766]	520
[524]	521	RETURN
	522	END
	523
	524
	525
	526
	527
	528
	529
	530
	531
	532
	533	c-----------------------------------------------------------------------------
	534	c Routine for reading SO4 data from files
	535	c-----------------------------------------------------------------------------
	536
	537
	538	SUBROUTINE getso4fromfile (cyr, so4)
	539	#include "netcdf.inc"
	540	#include "dimensions.h"
	541	#include "dimphy.h"
	542	CHARACTER*15 fname
	543	CHARACTER*4 cyr
	544
	545	CHARACTER*6 cvar
	546	INTEGER START(3), COUNT(3)
	547	INTEGER STATUS, NCID, VARID
	548	INTEGER imth, i, j, k, ny
	549	PARAMETER (ny=jjm+1)
	550
	551
[782]	552	REAL so4mth(iim, ny, klev)
	553	REAL so4(iim, ny, klev, 12)
[524]	554
	555
	556	fname = 'so4.run'//cyr//'.cdf'
	557
	558	write (,) 'reading ', fname
	559	STATUS = NF_OPEN (fname, NF_NOWRITE, NCID)
	560	IF (STATUS .NE. NF_NOERR) write (,) 'err in open ',status
	561
	562	DO imth=1, 12
	563	IF (imth.eq.1) THEN
	564	cvar='SO4JAN'
	565	ELSEIF (imth.eq.2) THEN
	566	cvar='SO4FEB'
	567	ELSEIF (imth.eq.3) THEN
	568	cvar='SO4MAR'
	569	ELSEIF (imth.eq.4) THEN
	570	cvar='SO4APR'
	571	ELSEIF (imth.eq.5) THEN
	572	cvar='SO4MAY'
	573	ELSEIF (imth.eq.6) THEN
	574	cvar='SO4JUN'
	575	ELSEIF (imth.eq.7) THEN
	576	cvar='SO4JUL'
	577	ELSEIF (imth.eq.8) THEN
	578	cvar='SO4AUG'
	579	ELSEIF (imth.eq.9) THEN
	580	cvar='SO4SEP'
	581	ELSEIF (imth.eq.10) THEN
	582	cvar='SO4OCT'
	583	ELSEIF (imth.eq.11) THEN
	584	cvar='SO4NOV'
	585	ELSEIF (imth.eq.12) THEN
	586	cvar='SO4DEC'
	587	ENDIF
	588	start(1)=1
	589	start(2)=1
	590	start(3)=1
	591	count(1)=iim
	592	count(2)=ny
	593	count(3)=klev
	594	c write(,) 'here i am'
	595	STATUS = NF_INQ_VARID (NCID, cvar, VARID)
	596	write (,) ncid,imth,cvar, varid
[782]	597
[524]	598	IF (STATUS .NE. NF_NOERR) write (,) 'err in read ',status
[782]	599
	600	#ifdef NC_DOUBLE
[801]	601	status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth)
[782]	602	#else
[801]	603	status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth)
[782]	604	#endif
[524]	605	IF (STATUS .NE. NF_NOERR) write (,) 'err in read data',status
	606
	607	DO k=1,klev
	608	DO j=1,jjm+1
	609	DO i=1,iim
	610	IF (so4mth(i,j,k).LT.0.) then
	611	write(,) 'this is shit'
	612	write(,) 'so4(',i,j,k,') =',so4mth(i,j,k)
	613	endif
	614	so4(i,j,k,imth)=so4mth(i,j,k)
	615	ENDDO
	616	ENDDO
	617	ENDDO
	618	ENDDO
	619
	620	STATUS = NF_CLOSE(NCID)
[782]	621	IF (STATUS .NE. NF_NOERR) write (,) 'err in closing file',status
	622
	623
[524]	624	END ! subroutine getso4fromfile
	625
	626
	627
	628
	629
	630
	631
	632
	633
	634
	635
	636
	637
	638
	639
	640

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

Download in other formats:

Original Format