/[lmdze]/trunk/phylmd/readsulfate_preind.f
ViewVC logotype

Diff of /trunk/phylmd/readsulfate_preind.f

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

revision 107 by guez, Thu Sep 4 10:40:24 2014 UTC revision 108 by guez, Tue Sep 16 14:00:41 2014 UTC
# Line 1  Line 1 
1  SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate)  module readsulfate_preind_m
2    
   ! Read in /calculate pre-industrial values of sulfate  
   
   use dimens_m  
   use dimphy  
   use temps  
   use SUPHEC_M  
   use chem  
   use getso4fromfile_m, only: getso4fromfile  
3    IMPLICIT none    IMPLICIT none
4    
5    ! Content:  contains
6    ! --------  
7    ! This routine reads in monthly mean values of sulfate aerosols and    SUBROUTINE readsulfate_preind(r_day, first, pi_sulfate)
8    ! returns a linearly interpolated daily-mean field.  
9    !      ! Read in /calculate pre-industrial values of sulfate
10    ! It does so for the preindustriel values of the sulfate, to a large part  
11    ! analogous to the routine readsulfate.      use dimens_m
12    !      use dimphy
13    ! Only Pb: Variables must be saved and don t have to be overwritten!      use temps
14    !      use SUPHEC_M
15    ! Author:      use chem
16    ! -------      use getso4fromfile_m, only: getso4fromfile
17    ! Johannes Quaas (quaas@lmd.jussieu.fr)  
18    ! 26/06/01      ! Content:
19    !      ! --------
20    ! Input:      ! This routine reads in monthly mean values of sulfate aerosols and
21    ! ------      ! returns a linearly interpolated daily-mean field.
22    double precision, intent(in)::  r_day                   ! Day of integration      !
23    LOGICAL, intent(in):: first                 ! First timestep      ! It does so for the preindustriel values of the sulfate, to a large part
24    ! (and therefore initialization necessary)      ! analogous to the routine readsulfate.
25    !      !
26    ! Output:      ! Only Pb: Variables must be saved and don t have to be overwritten!
27    ! -------      !
28    double precision  pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data,      ! Author:
29    !  from file)      ! -------
30    !      ! Johannes Quaas (quaas@lmd.jussieu.fr)
31    ! Local Variables:      ! 26/06/01
32    ! ----------------      !
33    INTEGER i, ig, k, it      ! Input:
34    INTEGER j, iday, ny, iyr      ! ------
35    parameter (ny=jjm+1)      real, intent(in)::  r_day                   ! Day of integration
36        LOGICAL, intent(in):: first                 ! First timestep
37    INTEGER im, day1, day2, im2      ! (and therefore initialization necessary)
38    double precision pi_so4_1(iim, jjm+1, klev, 12)      !
39        ! Output:
40    double precision pi_so4(klon, klev, 12)  ! SO4 in right dimension      ! -------
41    SAVE pi_so4      real pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data,
42    double precision pi_so4_out(klon, klev)      !  from file)
43    SAVE pi_so4_out      !
44        ! Local Variables:
45    CHARACTER(len=4) cyear      ! ----------------
46    LOGICAL lnewday      INTEGER i, ig, k, it
47        INTEGER j, iday, ny, iyr
48        parameter (ny=jjm+1)
49    
50    iday = INT(r_day)      INTEGER im, day1, day2, im2
51        double precision pi_so4_1(iim, jjm+1, klev, 12)
52    ! Get the year of the run  
53    iyr  = iday/360      double precision pi_so4(klon, klev, 12)  ! SO4 in right dimension
54        SAVE pi_so4
55    ! Get the day of the actual year:      double precision pi_so4_out(klon, klev)
56    iday = iday-iyr*360      SAVE pi_so4_out
57    
58    ! 0.02 is about 0.5/24, namly less than half an hour      CHARACTER(len=4) cyear
59    lnewday = (r_day-FLOAT(iday).LT.0.02)      LOGICAL lnewday
60    
61    ! ---------------------------------------------  
62    ! All has to be done only, if a new day begins!  
63    ! ---------------------------------------------      iday = INT(r_day)
64    
65    IF (lnewday.OR.first) THEN      ! Get the year of the run
66       im = iday/30 +1 ! the actual month      iyr  = iday/360
67    
68       ! annee_ref is the initial year (defined in temps.h)      ! Get the day of the actual year:
69       iyr = iyr + annee_ref      iday = iday-iyr*360
70    
71        ! 0.02 is about 0.5/24, namly less than half an hour
72       IF (first) THEN      lnewday = (r_day-FLOAT(iday).LT.0.02)
73          cyear='.nat'  
74          CALL getso4fromfile(cyear,pi_so4_1)      ! ---------------------------------------------
75        ! All has to be done only, if a new day begins!
76          ! Transform the horizontal 2D-field into the physics-field      ! ---------------------------------------------
77          ! (Also the levels and the latitudes have to be inversed)  
78        IF (lnewday.OR.first) THEN
79          ! Initialize field         im = iday/30 +1 ! the actual month
80          DO it=1,12  
81             DO k=1,klev         ! annee_ref is the initial year (defined in temps.h)
82                DO i=1,klon         iyr = iyr + annee_ref
83                   pi_so4(i,k,it)=0.  
84                ENDDO  
85             ENDDO         IF (first) THEN
86          ENDDO            cyear='.nat'
87              CALL getso4fromfile(cyear,pi_so4_1)
88          write (*,*) 'preind: finished reading', FLOAT(iim)  
89          DO it=1,12            ! Transform the horizontal 2D-field into the physics-field
90             DO k=1, klev            ! (Also the levels and the latitudes have to be inversed)
91                ! a) at the poles, use the zonal mean:  
92                DO i=1,iim            ! Initialize field
93                   ! North pole            DO it=1,12
94                   pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)               DO k=1,klev
95                   ! South pole                  DO i=1,klon
96                   pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)                     pi_so4(i,k,it)=0.
97                ENDDO                  ENDDO
98                pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)               ENDDO
99                pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)            ENDDO
100    
101                ! b) the values between the poles:            write (*,*) 'preind: finished reading', FLOAT(iim)
102                ig=1            DO it=1,12
103                DO j=2,jjm               DO k=1, klev
104                   DO i=1,iim                  ! a) at the poles, use the zonal mean:
105                      ig=ig+1                  DO i=1,iim
106                      if (ig.gt.klon) write (*,*) 'shit'                     ! North pole
107                      pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)                     pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
108                   ENDDO                     ! South pole
109                ENDDO                     pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
110                IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'                  ENDDO
111             ENDDO ! Loop over k (vertical)                  pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
112          ENDDO ! Loop over it (months)                  pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
113    
114       ENDIF                     ! Had to read new data?                  ! b) the values between the poles:
115                    ig=1
116                    DO j=2,jjm
117       ! Interpolate to actual day:                     DO i=1,iim
118       IF (iday.LT.im*30-15) THEN                        ig=ig+1
119          ! in the first half of the month use month before and actual month                        if (ig.gt.klon) write (*,*) 'shit'
120          im2=im-1                        pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)
121          day1 = im2*30+15                     ENDDO
122          day2 = im2*30-15                  ENDDO
123          IF (im2.LE.0) THEN                  IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
124             ! the month is january, thus the month before december               ENDDO ! Loop over k (vertical)
125             im2=12            ENDDO ! Loop over it (months)
126          ENDIF  
127          DO k=1,klev         ENDIF                     ! Had to read new data?
128             DO i=1,klon  
129                pi_sulfate(i,k) = pi_so4(i,k,im2)   &  
130                     - FLOAT(iday-day2)/FLOAT(day1-day2) &         ! Interpolate to actual day:
131                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))         IF (iday.LT.im*30-15) THEN
132                IF (pi_sulfate(i,k).LT.0.) THEN            ! in the first half of the month use month before and actual month
133                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2            im2=im-1
134                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &            day1 = im2*30+15
135                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &            day2 = im2*30-15
136                        pi_so4(i,k,im2) - pi_so4(i,k,im)            IF (im2.LE.0) THEN
137                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2               ! the month is january, thus the month before december
138                   stop 'pi_sulfate'               im2=12
139                endif            ENDIF
140             ENDDO            DO k=1,klev
141          ENDDO               DO i=1,klon
142       ELSE                  pi_sulfate(i,k) = pi_so4(i,k,im2)   &
143          ! the second half of the month                       - FLOAT(iday-day2)/FLOAT(day1-day2) &
144          im2=im+1                       * (pi_so4(i,k,im2) - pi_so4(i,k,im))
145          day1 = im*30+15                  IF (pi_sulfate(i,k).LT.0.) THEN
146          IF (im2.GT.12) THEN                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
147             ! the month is december, the following thus january                     IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
148             im2=1                          write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
149          ENDIF                          pi_so4(i,k,im2) - pi_so4(i,k,im)
150          day2 = im*30-15                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
151                       stop 'pi_sulfate'
152          DO k=1,klev                  endif
153             DO i=1,klon               ENDDO
154                pi_sulfate(i,k) = pi_so4(i,k,im2)   &            ENDDO
155                     - FLOAT(iday-day2)/FLOAT(day1-day2) &         ELSE
156                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))            ! the second half of the month
157                IF (pi_sulfate(i,k).LT.0.) THEN            im2=im+1
158                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2            day1 = im*30+15
159                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &            IF (im2.GT.12) THEN
160                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &               ! the month is december, the following thus january
161                        pi_so4(i,k,im2) - pi_so4(i,k,im)               im2=1
162                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2            ENDIF
163                   stop 'pi_sulfate'            day2 = im*30-15
164                endif  
165             ENDDO            DO k=1,klev
166          ENDDO               DO i=1,klon
167       ENDIF                  pi_sulfate(i,k) = pi_so4(i,k,im2)   &
168                         - FLOAT(iday-day2)/FLOAT(day1-day2) &
169                         * (pi_so4(i,k,im2) - pi_so4(i,k,im))
170       !JLD      ! The sulfate concentration [molec cm-3] is read in.                  IF (pi_sulfate(i,k).LT.0.) THEN
171       !JLD      ! Convert it into mass [ug SO4/m3]                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
172       !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]                     IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
173       DO k=1,klev                          write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
174          DO i=1,klon                          pi_so4(i,k,im2) - pi_so4(i,k,im)
175             !JLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
176             !JLD     .           /n_avogadro*1.e+12                     stop 'pi_sulfate'
177             pi_so4_out(i,k) = pi_sulfate(i,k)                  endif
178          ENDDO               ENDDO
179       ENDDO            ENDDO
180           ENDIF
181    ELSE ! If no new day, use old data:  
182       DO k=1,klev  
183          DO i=1,klon         !JLD      ! The sulfate concentration [molec cm-3] is read in.
184             pi_sulfate(i,k) = pi_so4_out(i,k)         !JLD      ! Convert it into mass [ug SO4/m3]
185          ENDDO         !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
186       ENDDO         DO k=1,klev
187    ENDIF ! Was this the beginning of a new day?            DO i=1,klon
188                 !JLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
189                 !JLD     .           /n_avogadro*1.e+12
190                 pi_so4_out(i,k) = pi_sulfate(i,k)
191              ENDDO
192           ENDDO
193    
194        ELSE ! If no new day, use old data:
195           DO k=1,klev
196              DO i=1,klon
197                 pi_sulfate(i,k) = pi_so4_out(i,k)
198              ENDDO
199           ENDDO
200        ENDIF ! Was this the beginning of a new day?
201    
202      END SUBROUTINE readsulfate_preind
203    
204  END SUBROUTINE readsulfate_preind  end module readsulfate_preind_m
Legend:
Removed from v.107  
changed lines
  Added in v.108
  ViewVC Help
Powered by ViewVC 1.1.21