source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/readsulfate.F @ 3919

Last change on this file since 3919 was 634, checked in by Laurent Fairhead, 20 years ago

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