Context Navigation

readsulfate.F @ 3722

Last change on this file since 3722 was 638, checked in by Laurent Fairhead, 19 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

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, iyr1, iyr2
54	parameter (ny=jjm+1)
55
56	INTEGER ismaller
57	CJLD 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
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:
117	! a) from actual 10-yr-period
118
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
128	ENDIF
129	WRITE(cyear,'(I4)') iyr1
130	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
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
138	WRITE(cyear,'(I4)') iyr2
139	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
140	CALL getso4fromfile(cyear, so4_2)
141
142	ENDIF
143
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)
150	. - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
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
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.
247	DO k=1,klev
248	DO i=1,klon
249	CJLD sulfate(i,k) = sulfate(i,k)*masse_so4
250	CJLD . /n_avogadro*1.e+12
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
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]
464	DO k=1,klev
465	DO i=1,klon
466	CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
467	CJLD . /n_avogadro*1.e+12
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