/[lmdze]/trunk/libf/phylmd/readsulfate.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/readsulfate.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
File size: 17304 byte(s)
Initial import
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/readsulfate.F,v 1.2 2005/05/19 08:27:15 fairhead Exp $
3     !
4     SUBROUTINE readsulfate (r_day, first, sulfate)
5    
6     use dimens_m
7     use dimphy
8     use temps
9     use YOMCST
10     IMPLICIT none
11    
12     c Content:
13     c --------
14     c This routine reads in monthly mean values of sulfate aerosols and
15     c returns a linearly interpolated dayly-mean field.
16     c
17     c
18     c Author:
19     c -------
20     c Johannes Quaas (quaas@lmd.jussieu.fr)
21     c 26/04/01
22     c
23     c Modifications:
24     c --------------
25     c 21/06/01: Make integrations of more than one year possible ;-)
26     c ATTENTION!! runs are supposed to start with Jan, 1. 1930
27     c (rday=1)
28     c
29     c 27/06/01: Correction: The model always has 360 days per year!
30     c 27/06/01: SO4 concentration rather than mixing ratio
31     c 27/06/01: 10yr-mean-values to interpolate
32     c 20/08/01: Correct the error through integer-values in interpolations
33     c 21/08/01: Introduce flag to read in just one decade
34     c
35     c Include-files:
36     c --------------
37     include "chem.h"
38     c
39     c Input:
40     c ------
41     REAL*8 r_day ! Day of integration
42     LOGICAL, intent(in):: 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     use dimens_m
282     use dimphy
283     use temps
284     use YOMCST
285     IMPLICIT none
286    
287     c Content:
288     c --------
289     c This routine reads in monthly mean values of sulfate aerosols and
290     c returns a linearly interpolated dayly-mean field.
291     c
292     c It does so for the preindustriel values of the sulfate, to a large part
293     c analogous to the routine readsulfate above.
294     c
295     c Only Pb: Variables must be saved and don t have to be overwritten!
296     c
297     c Author:
298     c -------
299     c Johannes Quaas (quaas@lmd.jussieu.fr)
300     c 26/06/01
301     c
302     c Modifications:
303     c --------------
304     c see above
305     c
306     c Include-files:
307     c --------------
308     include "chem.h"
309     c
310     c Input:
311     c ------
312     REAL*8 r_day ! Day of integration
313     LOGICAL, intent(in):: 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     use dimens_m
500     use dimphy
501     include "netcdf.inc"
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    

  ViewVC Help
Powered by ViewVC 1.1.21