source: LMDZ4/branches/V3_test/libf/phylmd/readsulfate.F @ 2238

Last change on this file since 2238 was 704, checked in by Laurent Fairhead, 18 years ago

Inclusion des modifs de Y. Meurdesoif pour la version V3
LF

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