source: LMDZ.3.3/branches/rel-LF/libf/phylmd/readsulfate.F @ 520

Last change on this file since 520 was 520, checked in by (none), 21 years ago

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