Context Navigation

readsulfate.F @ 5420

Last change on this file since 5420 was 638, checked in by Laurent Fairhead, 20 years ago
Synchro avec la version de LMDZ.3.3 pour Arnaud LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.8 KB

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

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

Download in other formats:

Original Format