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