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

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

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

trunk/phylmd/readsulfate_preind.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC trunk/Sources/phylmd/readsulfate_preind.f revision 188 by guez, Tue Mar 22 16:31:39 2016 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(dayvrai, time, 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 dynetat0_m, only: annee_ref
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)      integer, intent(in):: dayvrai
36        ! current day number, based at value 1 on January 1st of annee_ref
37    INTEGER im, day1, day2, im2, ismaller  
38    double precision pi_so4_1(iim, jjm+1, klev, 12)      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
39    
40    double precision pi_so4(klon, klev, 12)  ! SO4 in right dimension      LOGICAL, intent(in):: first                 ! First timestep
41    SAVE pi_so4      ! (and therefore initialization necessary)
42    double precision pi_so4_out(klon, klev)      !
43    SAVE pi_so4_out      ! Output:
44        ! -------
45    CHARACTER*4 cyear      real pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data,
46    LOGICAL lnewday      !  from file)
47        !
48        ! Local Variables:
49        ! ----------------
50    iday = INT(r_day)      INTEGER i, ig, k, it
51        INTEGER j, iday, iyr
52    ! Get the year of the run  
53    iyr  = iday/360      INTEGER im, day1, day2, im2
54        double precision pi_so4_1(iim, jjm+1, klev, 12)
55    ! Get the day of the actual year:  
56    iday = iday-iyr*360      double precision pi_so4(klon, klev, 12)  ! SO4 in right dimension
57        SAVE pi_so4
58    ! 0.02 is about 0.5/24, namly less than half an hour      double precision pi_so4_out(klon, klev)
59    lnewday = (r_day-FLOAT(iday).LT.0.02)      SAVE pi_so4_out
60    
61    ! ---------------------------------------------      CHARACTER(len=4) cyear
62    ! All has to be done only, if a new day begins!      LOGICAL lnewday
63    ! ---------------------------------------------  
64    
65    IF (lnewday.OR.first) THEN  
66       im = iday/30 +1 ! the actual month      iday = dayvrai
67    
68       ! annee_ref is the initial year (defined in temps.h)      ! Get the year of the run
69       iyr = iyr + annee_ref      iyr  = iday/360
70    
71        ! Get the day of the actual year:
72       IF (first) THEN      iday = iday-iyr*360
73          cyear='.nat'  
74          CALL getso4fromfile(cyear,pi_so4_1)      ! 0.02 is about 0.5/24, namly less than half an hour
75        lnewday = time < 0.02
76          ! Transform the horizontal 2D-field into the physics-field  
77          ! (Also the levels and the latitudes have to be inversed)      ! ---------------------------------------------
78        ! All has to be done only, if a new day begins!
79          ! Initialize field      ! ---------------------------------------------
80          DO it=1,12  
81             DO k=1,klev      IF (lnewday.OR.first) THEN
82                DO i=1,klon         im = iday/30 +1 ! the actual month
83                   pi_so4(i,k,it)=0.  
84                ENDDO         ! annee_ref is the initial year (defined in temps.h)
85             ENDDO         iyr = iyr + annee_ref
86          ENDDO  
87    
88          write (*,*) 'preind: finished reading', FLOAT(iim)         IF (first) THEN
89          DO it=1,12            cyear='.nat'
90             DO k=1, klev            CALL getso4fromfile(cyear,pi_so4_1)
91                ! a) at the poles, use the zonal mean:  
92                DO i=1,iim            ! Transform the horizontal 2D-field into the physics-field
93                   ! North pole            ! (Also the levels and the latitudes have to be inversed)
94                   pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)  
95                   ! South pole            ! Initialize field
96                   pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)            DO it=1,12
97                ENDDO               DO k=1,klev
98                pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)                  DO i=1,klon
99                pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)                     pi_so4(i,k,it)=0.
100                    ENDDO
101                ! b) the values between the poles:               ENDDO
102                ig=1            ENDDO
103                DO j=2,jjm  
104                   DO i=1,iim            write (*,*) 'preind: finished reading', FLOAT(iim)
105                      ig=ig+1            DO it=1,12
106                      if (ig.gt.klon) write (*,*) 'shit'               DO k=1, klev
107                      pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)                  ! a) at the poles, use the zonal mean:
108                   ENDDO                  DO i=1,iim
109                ENDDO                     ! North pole
110                IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'                     pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
111             ENDDO ! Loop over k (vertical)                     ! South pole
112          ENDDO ! Loop over it (months)                     pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
113                    ENDDO
114       ENDIF                     ! Had to read new data?                  pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
115                    pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
116    
117       ! Interpolate to actual day:                  ! b) the values between the poles:
118       IF (iday.LT.im*30-15) THEN                  ig=1
119          ! in the first half of the month use month before and actual month                  DO j=2,jjm
120          im2=im-1                     DO i=1,iim
121          day1 = im2*30+15                        ig=ig+1
122          day2 = im2*30-15                        if (ig.gt.klon) write (*,*) 'shit'
123          IF (im2.LE.0) THEN                        pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)
124             ! the month is january, thus the month before december                     ENDDO
125             im2=12                  ENDDO
126          ENDIF                  IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
127          DO k=1,klev               ENDDO ! Loop over k (vertical)
128             DO i=1,klon            ENDDO ! Loop over it (months)
129                pi_sulfate(i,k) = pi_so4(i,k,im2)   &  
130                     - FLOAT(iday-day2)/FLOAT(day1-day2) &         ENDIF                     ! Had to read new data?
131                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))  
132                IF (pi_sulfate(i,k).LT.0.) THEN  
133                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2         ! Interpolate to actual day:
134                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &         IF (iday.LT.im*30-15) THEN
135                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &            ! in the first half of the month use month before and actual month
136                        pi_so4(i,k,im2) - pi_so4(i,k,im)            im2=im-1
137                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2            day1 = im2*30+15
138                   stop 'pi_sulfate'            day2 = im2*30-15
139                endif            IF (im2.LE.0) THEN
140             ENDDO               ! the month is january, thus the month before december
141          ENDDO               im2=12
142       ELSE            ENDIF
143          ! the second half of the month            DO k=1,klev
144          im2=im+1               DO i=1,klon
145          day1 = im*30+15                  pi_sulfate(i,k) = pi_so4(i,k,im2)   &
146          IF (im2.GT.12) THEN                       - FLOAT(iday-day2)/FLOAT(day1-day2) &
147             ! the month is december, the following thus january                       * (pi_so4(i,k,im2) - pi_so4(i,k,im))
148             im2=1                  IF (pi_sulfate(i,k).LT.0.) THEN
149          ENDIF                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
150          day2 = im*30-15                     IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
151                            write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
152          DO k=1,klev                          pi_so4(i,k,im2) - pi_so4(i,k,im)
153             DO i=1,klon                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
154                pi_sulfate(i,k) = pi_so4(i,k,im2)   &                     stop 'pi_sulfate'
155                     - FLOAT(iday-day2)/FLOAT(day1-day2) &                  endif
156                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))               ENDDO
157                IF (pi_sulfate(i,k).LT.0.) THEN            ENDDO
158                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2         ELSE
159                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &            ! the second half of the month
160                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &            im2=im+1
161                        pi_so4(i,k,im2) - pi_so4(i,k,im)            day1 = im*30+15
162                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2            IF (im2.GT.12) THEN
163                   stop 'pi_sulfate'               ! the month is december, the following thus january
164                endif               im2=1
165             ENDDO            ENDIF
166          ENDDO            day2 = im*30-15
167       ENDIF  
168              DO k=1,klev
169                 DO i=1,klon
170       !JLD      ! The sulfate concentration [molec cm-3] is read in.                  pi_sulfate(i,k) = pi_so4(i,k,im2)   &
171       !JLD      ! Convert it into mass [ug SO4/m3]                       - FLOAT(iday-day2)/FLOAT(day1-day2) &
172       !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]                       * (pi_so4(i,k,im2) - pi_so4(i,k,im))
173       DO k=1,klev                  IF (pi_sulfate(i,k).LT.0.) THEN
174          DO i=1,klon                     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
175             !JLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4                     IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
176             !JLD     .           /n_avogadro*1.e+12                          write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
177             pi_so4_out(i,k) = pi_sulfate(i,k)                          pi_so4(i,k,im2) - pi_so4(i,k,im)
178          ENDDO                     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
179       ENDDO                     stop 'pi_sulfate'
180                    endif
181    ELSE ! If no new day, use old data:               ENDDO
182       DO k=1,klev            ENDDO
183          DO i=1,klon         ENDIF
184             pi_sulfate(i,k) = pi_so4_out(i,k)  
185          ENDDO  
186       ENDDO         !JLD      ! The sulfate concentration [molec cm-3] is read in.
187    ENDIF ! Was this the beginning of a new day?         !JLD      ! Convert it into mass [ug SO4/m3]
188           !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
189           DO k=1,klev
190              DO i=1,klon
191                 !JLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
192                 !JLD     .           /n_avogadro*1.e+12
193                 pi_so4_out(i,k) = pi_sulfate(i,k)
194              ENDDO
195           ENDDO
196    
197        ELSE ! If no new day, use old data:
198           DO k=1,klev
199              DO i=1,klon
200                 pi_sulfate(i,k) = pi_so4_out(i,k)
201              ENDDO
202           ENDDO
203        ENDIF ! Was this the beginning of a new day?
204    
205      END SUBROUTINE readsulfate_preind
206    
207  END SUBROUTINE readsulfate_preind  end module readsulfate_preind_m
Legend:
Removed from v.82  
changed lines
  Added in v.188
  ViewVC Help
Powered by ViewVC 1.1.21