/[lmdze]/trunk/Sources/phylmd/readsulfate.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/readsulfate.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 190 by guez, Tue Mar 22 16:31:39 2016 UTC revision 191 by guez, Mon May 9 19:56:28 2016 UTC
# Line 9  contains Line 9  contains
9      ! From LMDZ4/libf/phylmd/readsulfate.F, version 1.2 2005/05/19      ! From LMDZ4/libf/phylmd/readsulfate.F, version 1.2 2005/05/19
10      ! 08:27:15 fairhead      ! 08:27:15 fairhead
11    
12      ! This routine reads in monthly mean values of sulfate aerosols and      ! This routine reads monthly mean values of sulfate aerosols and
13      ! returns a linearly interpolated daily-mean field.      ! returns a linearly interpolated daily-mean field.
14    
15      ! Author: Johannes Quaas (quaas@lmd.jussieu.fr)      ! Author: Johannes Quaas (quaas@lmd.jussieu.fr)
16      ! 26/04/01      ! April 26th, 2001
17    
     ! Modifications:  
     ! 21/06/01: Make integrations of more than one year possible ;-)  
18      ! ATTENTION!! runs are supposed to start with Jan, 1. 1930      ! ATTENTION!! runs are supposed to start with Jan, 1. 1930
19      ! (rday=1)      ! (rday = 1)
20    
21      ! 27/06/01: Correction: The model always has 360 days per year!      ! The model always has 360 days per year.
22      ! 27/06/01: SO4 concentration rather than mixing ratio      ! SO4 concentration rather than mixing ratio
23      ! 27/06/01: 10yr-mean-values to interpolate      ! 10yr-mean-values to interpolate
24      ! 20/08/01: Correct the error through integer-values in interpolations      ! Introduce flag to read in just one decade
25      ! 21/08/01: Introduce flag to read in just one decade  
26        use dimens_m, only: iim, jjm
27      USE dimens_m, ONLY: iim, jjm      use dimphy, only: klon, klev
28      USE dimphy, ONLY: klev, klon      use dynetat0_m, only: annee_ref
     USE dynetat0_m, ONLY: annee_ref  
29      use getso4fromfile_m, only: getso4fromfile      use getso4fromfile_m, only: getso4fromfile
30    
     ! Input:  
   
31      integer, intent(in):: dayvrai      integer, intent(in):: dayvrai
32      ! current day number, based at value 1 on January 1st of annee_ref      ! current day number, based at value 1 on January 1st of annee_ref
33    
34      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
35    
36      LOGICAL, intent(in):: first                 ! First timestep      LOGICAL, intent(in):: first ! First timestep
37      ! (and therefore initialization necessary)      ! (and therefore initialization necessary)
38    
39      ! Output:      real, intent(out):: sulfate(klon, klev)
40        ! mass of sulfate (monthly mean data, from file) (micro g SO4 / m3)
     real sulfate (klon, klev)  ! Mass of sulfate (monthly mean data,  
     !  from file) [ug SO4/m3]  
   
     ! Local Variables:  
41    
42        ! Local:
43      INTEGER i, ig, k, it      INTEGER i, ig, k, it
44      INTEGER j, iday, iyr, iyr1, iyr2      INTEGER j, iday, iyr, iyr1, iyr2
45        CHARACTER(len = 4) cyear
     CHARACTER(len=4) cyear  
   
46      INTEGER im, day1, day2, im2      INTEGER im, day1, day2, im2
47      double precision so4_1(iim, jjm+1, klev, 12)      real so4_1(iim, jjm + 1, klev, 12)
48      double precision so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions      real so4_2(iim, jjm + 1, klev, 12) ! sulfate distributions
49        double precision, save:: so4(klon, klev, 12) ! SO4 in right dimension
50      double precision so4(klon, klev, 12)  ! SO4 in right dimension      double precision, save:: so4_out(klon, klev)
     SAVE so4  
     double precision so4_out(klon, klev)  
     SAVE so4_out  
   
51      LOGICAL lnewday      LOGICAL lnewday
     LOGICAL lonlyone  
     PARAMETER (lonlyone=.FALSE.)  
52    
53      !--------------------------------------------------------------------      !----------------------------------------------------------------
54    
55      iday = dayvrai      iday = dayvrai
56    
57      ! Get the year of the run      ! Get the year of the run
58      iyr  = iday/360      iyr = iday/360
59    
60      ! Get the day of the actual year:      ! Get the day of the actual year:
61      iday = iday-iyr*360      iday = iday - iyr*360
62    
63      ! 0.02 is about 0.5/24, namly less than half an hour      ! 0.02 is about 0.5/24, namly less than half an hour
64      lnewday = time < 0.02      lnewday = time < 0.02
65    
66      ! All has to be done only, if a new day begins!      ! All has to be done only, if a new day begins
67    
68        test_newday: IF (lnewday .OR. first) THEN
69           im = iday/30 + 1 ! the actual month
70    
     IF (lnewday.OR.first) THEN  
        im = iday/30 +1 ! the actual month  
71         ! annee_ref is the initial year (defined in temps.h)         ! annee_ref is the initial year (defined in temps.h)
72         iyr = iyr + annee_ref         iyr = iyr + annee_ref
73    
74         ! Do I have to read new data? (Is this the first day of a year?)         ! Do I have to read new data? (Is this the first day of a year?)
75         IF (first.OR.iday.EQ.1.) THEN         IF (first .OR. iday == 1.) THEN
76            ! Initialize values            ! Initialize field
77            DO it=1,12            DO it = 1, 12
78               DO k=1,klev               DO k = 1, klev
79                  DO i=1,klon                  DO i = 1, klon
80                     so4(i,k,it)=0.                     so4(i, k, it) = 0.
81                  ENDDO                  ENDDO
82               ENDDO               ENDDO
83            ENDDO            ENDDO
84    
85            IF (iyr .lt. 1850) THEN            IF (iyr < 1850) THEN
86               cyear='.nat'               cyear = '.nat'
87               WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear               print *, 'getso4 iyr = ', iyr, ' ', cyear
88               CALL getso4fromfile(cyear, so4_1)               CALL getso4fromfile(cyear, so4_1)
89            ELSE IF (iyr .ge. 2100) THEN            ELSE IF (iyr >= 2100) THEN
90               cyear='2100'               cyear = '2100'
91               WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear               print *, 'getso4 iyr = ', iyr, ' ', cyear
92               CALL getso4fromfile(cyear, so4_1)               CALL getso4fromfile(cyear, so4_1)
93            ELSE            ELSE
   
94               ! Read in data:               ! Read in data:
95               ! a) from actual 10-yr-period               ! a) from actual 10-yr-period
96    
97               IF (iyr.LT.1900) THEN               IF (iyr < 1900) THEN
98                  iyr1 = 1850                  iyr1 = 1850
99                  iyr2 = 1900                  iyr2 = 1900
100               ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN               ELSE IF (iyr >= 1900.and.iyr < 1920) THEN
101                  iyr1 = 1900                  iyr1 = 1900
102                  iyr2 = 1920                  iyr2 = 1920
103               ELSE               ELSE
104                  iyr1 = INT(iyr/10)*10                  iyr1 = INT(iyr/10)*10
105                  iyr2 = INT(1+iyr/10)*10                  iyr2 = INT(1 + iyr/10)*10
106               ENDIF               ENDIF
107               WRITE(cyear,'(I4)') iyr1               WRITE(cyear, '(I4)') iyr1
108               WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear               print *, 'getso4 iyr = ', iyr, ' ', cyear
109               CALL getso4fromfile(cyear, so4_1)               CALL getso4fromfile(cyear, so4_1)
110    
111               ! If to read two decades:               ! Read two decades:
112               IF (.NOT.lonlyone) THEN               ! b) from the next following one
113                 WRITE(cyear, '(I4)') iyr2
114                  ! b) from the next following one               print *, 'getso4 iyr = ', iyr, ' ', cyear
115                  WRITE(cyear,'(I4)') iyr2               CALL getso4fromfile(cyear, so4_2)
                 WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear  
                 CALL getso4fromfile(cyear, so4_2)  
   
              ENDIF  
116    
117               ! Interpolate linarily to the actual year:               ! Interpolate linarily to the actual year:
118               DO it=1,12               DO it = 1, 12
119                  DO k=1,klev                  DO k = 1, klev
120                     DO j=1,jjm                     DO j = 1, jjm
121                        DO i=1,iim                        DO i = 1, iim
122                           so4_1(i,j,k,it)=so4_1(i,j,k,it) &                           so4_1(i, j, k, it) = so4_1(i, j, k, it) &
123                                - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &                                - REAL(iyr - iyr1)/REAL(iyr2 - iyr1) &
124                                * (so4_1(i,j,k,it) - so4_2(i,j,k,it))                                * (so4_1(i, j, k, it) - so4_2(i, j, k, it))
125                        ENDDO                        ENDDO
126                     ENDDO                     ENDDO
127                  ENDDO                  ENDDO
128               ENDDO               ENDDO
129              ENDIF
           ENDIF !lonlyone  
130    
131            ! Transform the horizontal 2D-field into the physics-field            ! Transform the horizontal 2D-field into the physics-field
132            ! (Also the levels and the latitudes have to be inversed)            ! (Also the levels and the latitudes have to be inversed)
133    
134            DO it=1,12            DO it = 1, 12
135               DO k=1, klev               DO k = 1, klev
136                  ! a) at the poles, use the zonal mean:                  ! a) at the poles, use the zonal mean:
137                  DO i=1,iim                  DO i = 1, iim
138                     ! North pole                     ! North pole
139                     so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)                     so4(1, k, it) = so4(1, k, it) &
140                            + so4_1(i, jjm + 1, klev + 1 - k, it)
141                     ! South pole                     ! South pole
142                     so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)                     so4(klon, k, it) = so4(klon, k, it) &
143                            + so4_1(i, 1, klev + 1 - k, it)
144                  ENDDO                  ENDDO
145                  so4(1,k,it)=so4(1,k,it)/FLOAT(iim)                  so4(1, k, it) = so4(1, k, it)/REAL(iim)
146                  so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)                  so4(klon, k, it) = so4(klon, k, it)/REAL(iim)
147    
148                  ! b) the values between the poles:                  ! b) the values between the poles:
149                  ig=1                  ig = 1
150                  DO j=2,jjm                  DO j = 2, jjm
151                     DO i=1,iim                     DO i = 1, iim
152                        ig=ig+1                        ig = ig + 1
153                        if (ig.gt.klon) write (*,*) 'shit'                        if (ig > klon) stop 1
154                        so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)                        so4(ig, k, it) = so4_1(i, jjm + 1 - j, klev + 1 - k, it)
155                     ENDDO                     ENDDO
156                  ENDDO                  ENDDO
157                  IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'                  IF (ig /= klon - 1) then
158                       print *, 'Error in readsulfate (var conversion)'
159                       STOP 1
160                    end IF
161               ENDDO ! Loop over k (vertical)               ENDDO ! Loop over k (vertical)
162            ENDDO ! Loop over it (months)            ENDDO ! Loop over it (months)
   
163         ENDIF ! Had to read new data?         ENDIF ! Had to read new data?
164    
165         ! Interpolate to actual day:         ! Interpolate to actual day:
166         IF (iday.LT.im*30-15) THEN         IF (iday < im*30 - 15) THEN
167            ! in the first half of the month use month before and actual month            ! in the first half of the month use month before and actual month
168            im2=im-1            im2 = im - 1
169            day2 = im2*30-15            day1 = im2*30 + 15
170            day1 = im2*30+15            day2 = im2*30 - 15
171            IF (im2.LE.0) THEN            IF (im2 <= 0) THEN
172               ! the month is january, thus the month before december               ! the month is january, thus the month before december
173               im2=12               im2 = 12
174            ENDIF            ENDIF
175            DO k=1,klev            DO k = 1, klev
176               DO i=1,klon               DO i = 1, klon
177                  sulfate(i,k) = so4(i,k,im2)   &                  sulfate(i, k) = so4(i, k, im2) &
178                       - FLOAT(iday-day2)/FLOAT(day1-day2) &                       - REAL(iday - day2)/REAL(day1 - day2) &
179                       * (so4(i,k,im2) - so4(i,k,im))                       * (so4(i, k, im2) - so4(i, k, im))
180                  IF (sulfate(i,k).LT.0.) THEN                  IF (sulfate(i, k) < 0.) THEN
181                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2                     IF (iday - day2 < 0.) write(*, *) 'iday - day2', iday - day2
182                     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &                     IF (so4(i, k, im2) - so4(i, k, im) < 0.) &
183                          write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &                          write(*, *) 'so4(i, k, im2) - so4(i, k, im)', &
184                          so4(i,k,im2) - so4(i,k,im)                          so4(i, k, im2) - so4(i, k, im)
185                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2                     IF (day1 - day2 < 0.) write(*, *) 'day1 - day2', day1 - day2
186                     stop 'sulfate'                     stop 1
187                  endif                  endif
188               ENDDO               ENDDO
189            ENDDO            ENDDO
190         ELSE         ELSE
191            ! the second half of the month            ! the second half of the month
192            im2=im+1            im2 = im + 1
193            IF (im2.GT.12) THEN            IF (im2 > 12) THEN
194               ! the month is december, the following thus january               ! the month is december, the following thus january
195               im2=1               im2 = 1
196            ENDIF            ENDIF
197            day2 = im*30-15            day2 = im*30 - 15
198            day1 = im*30+15            day1 = im*30 + 15
199            DO k=1,klev            DO k = 1, klev
200               DO i=1,klon               DO i = 1, klon
201                  sulfate(i,k) = so4(i,k,im2)   &                  sulfate(i, k) = so4(i, k, im2) &
202                       - FLOAT(iday-day2)/FLOAT(day1-day2) &                       - REAL(iday - day2)/REAL(day1 - day2) &
203                       * (so4(i,k,im2) - so4(i,k,im))                       * (so4(i, k, im2) - so4(i, k, im))
204                  IF (sulfate(i,k).LT.0.) THEN                  IF (sulfate(i, k) < 0.) THEN
205                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2                     IF (iday - day2 < 0.) write(*, *) 'iday - day2', iday - day2
206                     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &                     IF (so4(i, k, im2) - so4(i, k, im) < 0.) &
207                          write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &                          write(*, *) 'so4(i, k, im2) - so4(i, k, im)', &
208                          so4(i,k,im2) - so4(i,k,im)                          so4(i, k, im2) - so4(i, k, im)
209                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2                     IF (day1 - day2 < 0.) write(*, *) 'day1 - day2', day1 - day2
210                     stop 'sulfate'                     stop 1
211                  endif                  endif
212               ENDDO               ENDDO
213            ENDDO            ENDDO
214         ENDIF         ENDIF
215    
216         !JLD      ! The sulfate concentration [molec cm-3] is read in.         DO k = 1, klev
217         !JLD      ! Convert it into mass [ug SO4/m3]            DO i = 1, klon
218         !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]               so4_out(i, k) = sulfate(i, k)
        ! The sulfate mass [ug SO4/m3] is read in.  
        DO k=1,klev  
           DO i=1,klon  
              !JLD            sulfate(i,k) = sulfate(i,k)*masse_so4  
              !JLD     .           /n_avogadro*1.e+12  
              so4_out(i,k) = sulfate(i,k)  
              IF (so4_out(i,k).LT.0)  &  
                   stop 'WAS SOLL DER SCHEISS ? '  
219            ENDDO            ENDDO
220         ENDDO         ENDDO
221      ELSE ! if no new day, use old data:      ELSE
222         DO k=1,klev         ! If no new day, use old data:
223            DO i=1,klon         DO k = 1, klev
224               sulfate(i,k) = so4_out(i,k)            DO i = 1, klon
225               IF (so4_out(i,k).LT.0)  &               sulfate(i, k) = so4_out(i, k)
                   stop 'WAS SOLL DER SCHEISS ? '  
226            ENDDO            ENDDO
227         ENDDO         ENDDO
228      ENDIF ! Did I have to do anything (was it a new day?)      ENDIF test_newday
229    
230    END SUBROUTINE readsulfate    END SUBROUTINE readsulfate
231    
Legend:
Removed from v.190  
changed lines
  Added in v.191
  ViewVC Help
Powered by ViewVC 1.1.21