/[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

trunk/libf/phylmd/readsulfate.f90 revision 68 by guez, Wed Nov 14 16:59:30 2012 UTC trunk/phylmd/readsulfate.f revision 129 by guez, Fri Feb 13 18:22:38 2015 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 dynetat0_m, ONLY: annee_ref
32    LOGICAL, intent(in):: first                 ! First timestep      use getso4fromfile_m, only: getso4fromfile
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)      CHARACTER(len=4) cyear
52    double precision so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions  
53        INTEGER im, day1, day2, im2
54    double precision so4(klon, klev, 12)  ! SO4 in right dimension      double precision so4_1(iim, jjm+1, klev, 12)
55    SAVE so4      double precision so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions
56    double precision so4_out(klon, klev)  
57    SAVE so4_out      double precision so4(klon, klev, 12)  ! SO4 in right dimension
58        SAVE so4
59    LOGICAL lnewday      double precision so4_out(klon, klev)
60    LOGICAL lonlyone      SAVE so4_out
61    PARAMETER (lonlyone=.FALSE.)  
62        LOGICAL lnewday
63    iday = INT(r_day)      LOGICAL lonlyone
64        PARAMETER (lonlyone=.FALSE.)
65    ! Get the year of the run  
66    iyr  = iday/360      !--------------------------------------------------------------------
67    
68    ! Get the day of the actual year:      iday = INT(r_day)
69    iday = iday-iyr*360  
70        ! Get the year of the run
71    ! 0.02 is about 0.5/24, namly less than half an hour      iyr  = iday/360
72    lnewday = (r_day-FLOAT(iday).LT.0.02)  
73        ! Get the day of the actual year:
74    ! ---------------------------------------------      iday = iday-iyr*360
75    ! All has to be done only, if a new day begins!  
76    ! ---------------------------------------------      ! 0.02 is about 0.5/24, namly less than half an hour
77        lnewday = (r_day-FLOAT(iday).LT.0.02)
78    IF (lnewday.OR.first) THEN  
79       im = iday/30 +1 ! the actual month      ! All has to be done only, if a new day begins!
80       ! annee_ref is the initial year (defined in temps.h)  
81       iyr = iyr + annee_ref      IF (lnewday.OR.first) THEN
82           im = iday/30 +1 ! the actual month
83       ! Do I have to read new data? (Is this the first day of a year?)         ! annee_ref is the initial year (defined in temps.h)
84       IF (first.OR.iday.EQ.1.) THEN         iyr = iyr + annee_ref
85          ! Initialize values  
86          DO it=1,12         ! Do I have to read new data? (Is this the first day of a year?)
87             DO k=1,klev         IF (first.OR.iday.EQ.1.) THEN
88                DO i=1,klon            ! Initialize values
89                   so4(i,k,it)=0.            DO it=1,12
90                ENDDO               DO k=1,klev
91             ENDDO                  DO i=1,klon
92          ENDDO                     so4(i,k,it)=0.
93                    ENDDO
94                 ENDDO
95          IF (iyr .lt. 1850) THEN            ENDDO
96             cyear='.nat'  
97             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear            IF (iyr .lt. 1850) THEN
98             CALL getso4fromfile(cyear, so4_1)               cyear='.nat'
99          ELSE IF (iyr .ge. 2100) THEN               WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
100             cyear='2100'               CALL getso4fromfile(cyear, so4_1)
101             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear            ELSE IF (iyr .ge. 2100) THEN
102             CALL getso4fromfile(cyear, so4_1)               cyear='2100'
103          ELSE               WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
104                 CALL getso4fromfile(cyear, so4_1)
105             ! Read in data:            ELSE
106             ! a) from actual 10-yr-period  
107                 ! Read in data:
108             IF (iyr.LT.1900) THEN               ! a) from actual 10-yr-period
109                iyr1 = 1850  
110                iyr2 = 1900               IF (iyr.LT.1900) THEN
111             ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN                  iyr1 = 1850
112                iyr1 = 1900                  iyr2 = 1900
113                iyr2 = 1920               ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
114             ELSE                  iyr1 = 1900
115                iyr1 = INT(iyr/10)*10                  iyr2 = 1920
116                iyr2 = INT(1+iyr/10)*10               ELSE
117             ENDIF                  iyr1 = INT(iyr/10)*10
118             WRITE(cyear,'(I4)') iyr1                  iyr2 = INT(1+iyr/10)*10
119             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear               ENDIF
120             CALL getso4fromfile(cyear, so4_1)               WRITE(cyear,'(I4)') iyr1
121                 WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
122                 CALL getso4fromfile(cyear, so4_1)
123             ! If to read two decades:  
124             IF (.NOT.lonlyone) THEN               ! If to read two decades:
125                 IF (.NOT.lonlyone) THEN
126                ! b) from the next following one  
127                WRITE(cyear,'(I4)') iyr2                  ! b) from the next following one
128                WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear                  WRITE(cyear,'(I4)') iyr2
129                CALL getso4fromfile(cyear, so4_2)                  WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
130                    CALL getso4fromfile(cyear, so4_2)
131             ENDIF  
132                 ENDIF
133             ! Interpolate linarily to the actual year:  
134             DO it=1,12               ! Interpolate linarily to the actual year:
135                DO k=1,klev               DO it=1,12
136                   DO j=1,jjm                  DO k=1,klev
137                      DO i=1,iim                     DO j=1,jjm
138                         so4_1(i,j,k,it)=so4_1(i,j,k,it) &                        DO i=1,iim
139                              - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &                           so4_1(i,j,k,it)=so4_1(i,j,k,it) &
140                              * (so4_1(i,j,k,it) - so4_2(i,j,k,it))                                - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &
141                      ENDDO                                * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
142                   ENDDO                        ENDDO
143                ENDDO                     ENDDO
144             ENDDO                  ENDDO
145                 ENDDO
146          ENDIF !lonlyone  
147              ENDIF !lonlyone
148          ! Transform the horizontal 2D-field into the physics-field  
149          ! (Also the levels and the latitudes have to be inversed)            ! Transform the horizontal 2D-field into the physics-field
150              ! (Also the levels and the latitudes have to be inversed)
151          DO it=1,12  
152             DO k=1, klev            DO it=1,12
153                ! a) at the poles, use the zonal mean:               DO k=1, klev
154                DO i=1,iim                  ! a) at the poles, use the zonal mean:
155                   ! North pole                  DO i=1,iim
156                   so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)                     ! North pole
157                   ! South pole                     so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
158                   so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)                     ! South pole
159                ENDDO                     so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
160                so4(1,k,it)=so4(1,k,it)/FLOAT(iim)                  ENDDO
161                so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)                  so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
162                    so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
163                ! b) the values between the poles:  
164                ig=1                  ! b) the values between the poles:
165                DO j=2,jjm                  ig=1
166                   DO i=1,iim                  DO j=2,jjm
167                      ig=ig+1                     DO i=1,iim
168                      if (ig.gt.klon) write (*,*) 'shit'                        ig=ig+1
169                      so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)                        if (ig.gt.klon) write (*,*) 'shit'
170                   ENDDO                        so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
171                ENDDO                     ENDDO
172                IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'                  ENDDO
173             ENDDO ! Loop over k (vertical)                  IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
174          ENDDO ! Loop over it (months)               ENDDO ! Loop over k (vertical)
175              ENDDO ! Loop over it (months)
176    
177       ENDIF ! Had to read new data?         ENDIF ! Had to read new data?
178    
179           ! Interpolate to actual day:
180       ! Interpolate to actual day:         IF (iday.LT.im*30-15) THEN
181       IF (iday.LT.im*30-15) THEN            ! in the first half of the month use month before and actual month
182          ! in the first half of the month use month before and actual month            im2=im-1
183          im2=im-1            day2 = im2*30-15
184          day2 = im2*30-15            day1 = im2*30+15
185          day1 = im2*30+15            IF (im2.LE.0) THEN
186          IF (im2.LE.0) THEN               ! the month is january, thus the month before december
187             ! the month is january, thus the month before december               im2=12
188             im2=12            ENDIF
189          ENDIF            DO k=1,klev
190          DO k=1,klev               DO i=1,klon
191             DO i=1,klon                  sulfate(i,k) = so4(i,k,im2)   &
192                sulfate(i,k) = so4(i,k,im2)   &                       - FLOAT(iday-day2)/FLOAT(day1-day2) &
193                     - FLOAT(iday-day2)/FLOAT(day1-day2) &                       * (so4(i,k,im2) - so4(i,k,im))
194                     * (so4(i,k,im2) - so4(i,k,im))                  IF (sulfate(i,k).LT.0.) THEN
195                IF (sulfate(i,k).LT.0.) THEN                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
196                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2                     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
197                   IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &                          write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
198                        write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &                          so4(i,k,im2) - so4(i,k,im)
199                        so4(i,k,im2) - so4(i,k,im)                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
200                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2                     stop 'sulfate'
201                   stop 'sulfate'                  endif
202                endif               ENDDO
203             ENDDO            ENDDO
204          ENDDO         ELSE
205       ELSE            ! the second half of the month
206          ! the second half of the month            im2=im+1
207          im2=im+1            IF (im2.GT.12) THEN
208          IF (im2.GT.12) THEN               ! the month is december, the following thus january
209             ! the month is december, the following thus january               im2=1
210             im2=1            ENDIF
211          ENDIF            day2 = im*30-15
212          day2 = im*30-15            day1 = im*30+15
213          day1 = im*30+15            DO k=1,klev
214          DO k=1,klev               DO i=1,klon
215             DO i=1,klon                  sulfate(i,k) = so4(i,k,im2)   &
216                sulfate(i,k) = so4(i,k,im2)   &                       - FLOAT(iday-day2)/FLOAT(day1-day2) &
217                     - FLOAT(iday-day2)/FLOAT(day1-day2) &                       * (so4(i,k,im2) - so4(i,k,im))
218                     * (so4(i,k,im2) - so4(i,k,im))                  IF (sulfate(i,k).LT.0.) THEN
219                IF (sulfate(i,k).LT.0.) THEN                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
220                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2                     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
221                   IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &                          write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
222                        write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &                          so4(i,k,im2) - so4(i,k,im)
223                        so4(i,k,im2) - so4(i,k,im)                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
224                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2                     stop 'sulfate'
225                   stop 'sulfate'                  endif
226                endif               ENDDO
227             ENDDO            ENDDO
228          ENDDO         ENDIF
229       ENDIF  
230           !JLD      ! The sulfate concentration [molec cm-3] is read in.
231           !JLD      ! Convert it into mass [ug SO4/m3]
232       !JLD      ! The sulfate concentration [molec cm-3] is read in.         !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
233       !JLD      ! Convert it into mass [ug SO4/m3]         ! The sulfate mass [ug SO4/m3] is read in.
234       !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]         DO k=1,klev
235       ! The sulfate mass [ug SO4/m3] is read in.            DO i=1,klon
236       DO k=1,klev               !JLD            sulfate(i,k) = sulfate(i,k)*masse_so4
237          DO i=1,klon               !JLD     .           /n_avogadro*1.e+12
238             !JLD            sulfate(i,k) = sulfate(i,k)*masse_so4               so4_out(i,k) = sulfate(i,k)
239             !JLD     .           /n_avogadro*1.e+12               IF (so4_out(i,k).LT.0)  &
240             so4_out(i,k) = sulfate(i,k)                    stop 'WAS SOLL DER SCHEISS ? '
241             IF (so4_out(i,k).LT.0)  &            ENDDO
242                  stop 'WAS SOLL DER SCHEISS ? '         ENDDO
243          ENDDO      ELSE ! if no new day, use old data:
244       ENDDO         DO k=1,klev
245    ELSE ! if no new day, use old data:            DO i=1,klon
246       DO k=1,klev               sulfate(i,k) = so4_out(i,k)
247          DO i=1,klon               IF (so4_out(i,k).LT.0)  &
248             sulfate(i,k) = so4_out(i,k)                    stop 'WAS SOLL DER SCHEISS ? '
249             IF (so4_out(i,k).LT.0)  &            ENDDO
250                  stop 'WAS SOLL DER SCHEISS ? '         ENDDO
251          ENDDO      ENDIF ! Did I have to do anything (was it a new day?)
252       ENDDO  
253    ENDIF ! Did I have to do anything (was it a new day?)    END SUBROUTINE readsulfate
254    
255  END SUBROUTINE readsulfate  end module readsulfate_m
Legend:
Removed from v.68  
changed lines
  Added in v.129
  ViewVC Help
Powered by ViewVC 1.1.21