source: LMDZ4/trunk/libf/phy_IPCC_AR4/readsulfate.F @ 1071

Last change on this file since 1071 was 868, checked in by Laurent Fairhead, 17 years ago

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