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