Context Navigation

readsulfate.F @ 558

Last change on this file since 558 was 523, checked in by lmdzadmin, 21 years ago

Dernieres modifs (ha, ha) sur le tag IPSL-CM4_IPCC, MED & JLD:

avant 1930 on force la lecture des sulfates naturels
bug de debordement de tableau dans physiq.F
taucalv passe à 10 ans

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

File size: 16.5 KB

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

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

Download in other formats:

Original Format