Context Navigation

readsulfate.F @ 4667

Last change on this file since 4667 was 524, checked in by lmdzadmin, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.3 KB

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

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

Download in other formats:

Original Format