source: LMDZ4/branches/LMDZ4_V3_patches/libf/phylmd/readsulfate.F @ 5404

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

Pb de lecture de la variable AC
LF

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