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