trunk/phylmd/readsulfate.f

module readsulfate_m

  IMPLICIT none

contains

  SUBROUTINE readsulfate(r_day, first, sulfate)

    ! From LMDZ4/libf/phylmd/readsulfate.F, version 1.2 2005/05/19
    ! 08:27:15 fairhead

    ! This routine reads in monthly mean values of sulfate aerosols and
    ! returns a linearly interpolated daily-mean field.

    ! Author: Johannes Quaas (quaas@lmd.jussieu.fr)
    ! 26/04/01

    ! Modifications:
    ! 21/06/01: Make integrations of more than one year possible ;-)
    ! ATTENTION!! runs are supposed to start with Jan, 1. 1930
    ! (rday=1)

    ! 27/06/01: Correction: The model always has 360 days per year!
    ! 27/06/01: SO4 concentration rather than mixing ratio
    ! 27/06/01: 10yr-mean-values to interpolate
    ! 20/08/01: Correct the error through integer-values in interpolations
    ! 21/08/01: Introduce flag to read in just one decade

    USE dimens_m, ONLY: iim, jjm
    USE dimphy, ONLY: klev, klon
    use getso4fromfile_m, only: getso4fromfile
    USE temps, ONLY: annee_ref

    ! Input:

    real, intent(in):: r_day                   ! Day of integration
    LOGICAL, intent(in):: first                 ! First timestep
    ! (and therefore initialization necessary)

    ! Output:

    real sulfate (klon, klev)  ! Mass of sulfate (monthly mean data,
    !  from file) [ug SO4/m3]

    ! Local Variables:

    INTEGER i, ig, k, it
    INTEGER j, iday, ny, iyr, iyr1, iyr2
    parameter (ny=jjm+1)

    INTEGER ismaller
    !JLD      INTEGER idec1, idec2 ! The two decadal data read ini
    CHARACTER*4 cyear

    INTEGER im, day1, day2, im2
    double precision so4_1(iim, jjm+1, klev, 12)
    double precision so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions

    double precision so4(klon, klev, 12)  ! SO4 in right dimension
    SAVE so4
    double precision so4_out(klon, klev)
    SAVE so4_out

    LOGICAL lnewday
    LOGICAL lonlyone
    PARAMETER (lonlyone=.FALSE.)

    !--------------------------------------------------------------------

    iday = INT(r_day)

    ! Get the year of the run
    iyr  = iday/360

    ! Get the day of the actual year:
    iday = iday-iyr*360

    ! 0.02 is about 0.5/24, namly less than half an hour
    lnewday = (r_day-FLOAT(iday).LT.0.02)

    ! All has to be done only, if a new day begins!

    IF (lnewday.OR.first) THEN
       im = iday/30 +1 ! the actual month
       ! annee_ref is the initial year (defined in temps.h)
       iyr = iyr + annee_ref

       ! Do I have to read new data? (Is this the first day of a year?)
       IF (first.OR.iday.EQ.1.) THEN
          ! Initialize values
          DO it=1,12
             DO k=1,klev
                DO i=1,klon
                   so4(i,k,it)=0.
                ENDDO
             ENDDO
          ENDDO

          IF (iyr .lt. 1850) THEN
             cyear='.nat'
             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
             CALL getso4fromfile(cyear, so4_1)
          ELSE IF (iyr .ge. 2100) THEN
             cyear='2100'
             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
             CALL getso4fromfile(cyear, so4_1)
          ELSE

             ! Read in data:
             ! a) from actual 10-yr-period

             IF (iyr.LT.1900) THEN
                iyr1 = 1850
                iyr2 = 1900
             ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
                iyr1 = 1900
                iyr2 = 1920
             ELSE
                iyr1 = INT(iyr/10)*10
                iyr2 = INT(1+iyr/10)*10
             ENDIF
             WRITE(cyear,'(I4)') iyr1
             WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
             CALL getso4fromfile(cyear, so4_1)

             ! If to read two decades:
             IF (.NOT.lonlyone) THEN

                ! b) from the next following one
                WRITE(cyear,'(I4)') iyr2
                WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
                CALL getso4fromfile(cyear, so4_2)

             ENDIF

             ! Interpolate linarily to the actual year:
             DO it=1,12
                DO k=1,klev
                   DO j=1,jjm
                      DO i=1,iim
                         so4_1(i,j,k,it)=so4_1(i,j,k,it) &
                              - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &
                              * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO

          ENDIF !lonlyone

          ! Transform the horizontal 2D-field into the physics-field
          ! (Also the levels and the latitudes have to be inversed)

          DO it=1,12
             DO k=1, klev
                ! a) at the poles, use the zonal mean:
                DO i=1,iim
                   ! North pole
                   so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
                   ! South pole
                   so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
                ENDDO
                so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
                so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)

                ! b) the values between the poles:
                ig=1
                DO j=2,jjm
                   DO i=1,iim
                      ig=ig+1
                      if (ig.gt.klon) write (*,*) 'shit'
                      so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
                   ENDDO
                ENDDO
                IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
             ENDDO ! Loop over k (vertical)
          ENDDO ! Loop over it (months)

       ENDIF ! Had to read new data?

       ! Interpolate to actual day:
       IF (iday.LT.im*30-15) THEN
          ! in the first half of the month use month before and actual month
          im2=im-1
          day2 = im2*30-15
          day1 = im2*30+15
          IF (im2.LE.0) THEN
             ! the month is january, thus the month before december
             im2=12
          ENDIF
          DO k=1,klev
             DO i=1,klon
                sulfate(i,k) = so4(i,k,im2)   &
                     - FLOAT(iday-day2)/FLOAT(day1-day2) &
                     * (so4(i,k,im2) - so4(i,k,im))
                IF (sulfate(i,k).LT.0.) THEN
                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
                   IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
                        write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
                        so4(i,k,im2) - so4(i,k,im)
                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
                   stop 'sulfate'
                endif
             ENDDO
          ENDDO
       ELSE
          ! the second half of the month
          im2=im+1
          IF (im2.GT.12) THEN
             ! the month is december, the following thus january
             im2=1
          ENDIF
          day2 = im*30-15
          day1 = im*30+15
          DO k=1,klev
             DO i=1,klon
                sulfate(i,k) = so4(i,k,im2)   &
                     - FLOAT(iday-day2)/FLOAT(day1-day2) &
                     * (so4(i,k,im2) - so4(i,k,im))
                IF (sulfate(i,k).LT.0.) THEN
                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
                   IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
                        write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
                        so4(i,k,im2) - so4(i,k,im)
                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
                   stop 'sulfate'
                endif
             ENDDO
          ENDDO
       ENDIF

       !JLD      ! The sulfate concentration [molec cm-3] is read in.
       !JLD      ! Convert it into mass [ug SO4/m3]
       !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
       ! The sulfate mass [ug SO4/m3] is read in.
       DO k=1,klev
          DO i=1,klon
             !JLD            sulfate(i,k) = sulfate(i,k)*masse_so4
             !JLD     .           /n_avogadro*1.e+12
             so4_out(i,k) = sulfate(i,k)
             IF (so4_out(i,k).LT.0)  &
                  stop 'WAS SOLL DER SCHEISS ? '
          ENDDO
       ENDDO
    ELSE ! if no new day, use old data:
       DO k=1,klev
          DO i=1,klon
             sulfate(i,k) = so4_out(i,k)
             IF (so4_out(i,k).LT.0)  &
                  stop 'WAS SOLL DER SCHEISS ? '
          ENDDO
       ENDDO
    ENDIF ! Did I have to do anything (was it a new day?)

  END SUBROUTINE readsulfate

end module readsulfate_m
1	guez	69	module readsulfate_m
2	guez	3
3	guez	68	IMPLICIT none
4	guez	3
5	guez	69	contains
6	guez	3
7	guez	69	SUBROUTINE readsulfate(r_day, first, sulfate)
8	guez	3
9	guez	69	! From LMDZ4/libf/phylmd/readsulfate.F, version 1.2 2005/05/19
10			! 08:27:15 fairhead
11	guez	3
12	guez	69	! This routine reads in monthly mean values of sulfate aerosols and
13			! returns a linearly interpolated daily-mean field.
14	guez	3
15	guez	69	! Author: Johannes Quaas (quaas@lmd.jussieu.fr)
16			! 26/04/01
17	guez	3
18	guez	69	! Modifications:
19			! 21/06/01: Make integrations of more than one year possible ;-)
20			! ATTENTION!! runs are supposed to start with Jan, 1. 1930
21			! (rday=1)
22	guez	3
23	guez	69	! 27/06/01: Correction: The model always has 360 days per year!
24			! 27/06/01: SO4 concentration rather than mixing ratio
25			! 27/06/01: 10yr-mean-values to interpolate
26			! 20/08/01: Correct the error through integer-values in interpolations
27			! 21/08/01: Introduce flag to read in just one decade
28	guez	3
29	guez	69	USE dimens_m, ONLY: iim, jjm
30			USE dimphy, ONLY: klev, klon
31			use getso4fromfile_m, only: getso4fromfile
32			USE temps, ONLY: annee_ref
33	guez	3
34	guez	69	! Input:
35	guez	3
36	guez	69	real, intent(in):: r_day ! Day of integration
37			LOGICAL, intent(in):: first ! First timestep
38			! (and therefore initialization necessary)
39	guez	3
40	guez	69	! Output:
41	guez	3
42	guez	69	real sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
43			! from file) [ug SO4/m3]
44	guez	3
45	guez	69	! Local Variables:
46	guez	3
47	guez	69	INTEGER i, ig, k, it
48			INTEGER j, iday, ny, iyr, iyr1, iyr2
49			parameter (ny=jjm+1)
50	guez	3
51	guez	69	INTEGER ismaller
52			!JLD INTEGER idec1, idec2 ! The two decadal data read ini
53			CHARACTER*4 cyear
54	guez	3
55	guez	69	INTEGER im, day1, day2, im2
56			double precision so4_1(iim, jjm+1, klev, 12)
57			double precision so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions
58	guez	3
59	guez	69	double precision so4(klon, klev, 12) ! SO4 in right dimension
60			SAVE so4
61			double precision so4_out(klon, klev)
62			SAVE so4_out
63	guez	3
64	guez	69	LOGICAL lnewday
65			LOGICAL lonlyone
66			PARAMETER (lonlyone=.FALSE.)
67	guez	3
68	guez	69	!--------------------------------------------------------------------
69	guez	3
70	guez	69	iday = INT(r_day)
71	guez	3
72	guez	69	! Get the year of the run
73			iyr = iday/360
74	guez	3
75	guez	69	! Get the day of the actual year:
76			iday = iday-iyr*360
77	guez	3
78	guez	69	! 0.02 is about 0.5/24, namly less than half an hour
79			lnewday = (r_day-FLOAT(iday).LT.0.02)
80	guez	3
81	guez	69	! All has to be done only, if a new day begins!
82	guez	3
83	guez	69	IF (lnewday.OR.first) THEN
84			im = iday/30 +1 ! the actual month
85			! annee_ref is the initial year (defined in temps.h)
86			iyr = iyr + annee_ref
87	guez	3
88	guez	69	! Do I have to read new data? (Is this the first day of a year?)
89			IF (first.OR.iday.EQ.1.) THEN
90			! Initialize values
91			DO it=1,12
92			DO k=1,klev
93			DO i=1,klon
94			so4(i,k,it)=0.
95			ENDDO
96			ENDDO
97			ENDDO
98	guez	3
99	guez	69	IF (iyr .lt. 1850) THEN
100			cyear='.nat'
101			WRITE(,) 'getso4 iyr=', iyr,' ',cyear
102			CALL getso4fromfile(cyear, so4_1)
103			ELSE IF (iyr .ge. 2100) THEN
104			cyear='2100'
105			WRITE(,) 'getso4 iyr=', iyr,' ',cyear
106			CALL getso4fromfile(cyear, so4_1)
107			ELSE
108	guez	3
109	guez	69	! Read in data:
110			! a) from actual 10-yr-period
111	guez	3
112	guez	69	IF (iyr.LT.1900) THEN
113			iyr1 = 1850
114			iyr2 = 1900
115			ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
116			iyr1 = 1900
117			iyr2 = 1920
118			ELSE
119			iyr1 = INT(iyr/10)*10
120			iyr2 = INT(1+iyr/10)*10
121			ENDIF
122			WRITE(cyear,'(I4)') iyr1
123			WRITE(,) 'getso4 iyr=', iyr,' ',cyear
124			CALL getso4fromfile(cyear, so4_1)
125	guez	3
126	guez	69	! If to read two decades:
127			IF (.NOT.lonlyone) THEN
128	guez	3
129	guez	69	! b) from the next following one
130			WRITE(cyear,'(I4)') iyr2
131			WRITE(,) 'getso4 iyr=', iyr,' ',cyear
132			CALL getso4fromfile(cyear, so4_2)
133	guez	3
134	guez	69	ENDIF
135
136			! Interpolate linarily to the actual year:
137			DO it=1,12
138			DO k=1,klev
139			DO j=1,jjm
140			DO i=1,iim
141			so4_1(i,j,k,it)=so4_1(i,j,k,it) &
142			- FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &
143			* (so4_1(i,j,k,it) - so4_2(i,j,k,it))
144			ENDDO
145			ENDDO
146			ENDDO
147			ENDDO
148
149			ENDIF !lonlyone
150
151			! Transform the horizontal 2D-field into the physics-field
152			! (Also the levels and the latitudes have to be inversed)
153
154			DO it=1,12
155			DO k=1, klev
156			! a) at the poles, use the zonal mean:
157			DO i=1,iim
158			! North pole
159			so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
160			! South pole
161			so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
162			ENDDO
163			so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
164			so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
165
166			! b) the values between the poles:
167			ig=1
168			DO j=2,jjm
169			DO i=1,iim
170			ig=ig+1
171			if (ig.gt.klon) write (,) 'shit'
172			so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
173			ENDDO
174			ENDDO
175			IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
176			ENDDO ! Loop over k (vertical)
177			ENDDO ! Loop over it (months)
178
179			ENDIF ! Had to read new data?
180
181			! Interpolate to actual day:
182			IF (iday.LT.im*30-15) THEN
183			! in the first half of the month use month before and actual month
184			im2=im-1
185			day2 = im2*30-15
186			day1 = im2*30+15
187			IF (im2.LE.0) THEN
188			! the month is january, thus the month before december
189			im2=12
190			ENDIF
191			DO k=1,klev
192			DO i=1,klon
193			sulfate(i,k) = so4(i,k,im2) &
194			- FLOAT(iday-day2)/FLOAT(day1-day2) &
195			* (so4(i,k,im2) - so4(i,k,im))
196			IF (sulfate(i,k).LT.0.) THEN
197			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
198			IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
199			write(,) 'so4(i,k,im2) - so4(i,k,im)', &
200			so4(i,k,im2) - so4(i,k,im)
201			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
202			stop 'sulfate'
203			endif
204			ENDDO
205			ENDDO
206			ELSE
207			! the second half of the month
208			im2=im+1
209			IF (im2.GT.12) THEN
210			! the month is december, the following thus january
211			im2=1
212			ENDIF
213			day2 = im*30-15
214			day1 = im*30+15
215			DO k=1,klev
216			DO i=1,klon
217			sulfate(i,k) = so4(i,k,im2) &
218			- FLOAT(iday-day2)/FLOAT(day1-day2) &
219			* (so4(i,k,im2) - so4(i,k,im))
220			IF (sulfate(i,k).LT.0.) THEN
221			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
222			IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
223			write(,) 'so4(i,k,im2) - so4(i,k,im)', &
224			so4(i,k,im2) - so4(i,k,im)
225			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
226			stop 'sulfate'
227			endif
228			ENDDO
229			ENDDO
230			ENDIF
231
232			!JLD ! The sulfate concentration [molec cm-3] is read in.
233			!JLD ! Convert it into mass [ug SO4/m3]
234			!JLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
235			! The sulfate mass [ug SO4/m3] is read in.
236			DO k=1,klev
237			DO i=1,klon
238			!JLD sulfate(i,k) = sulfate(i,k)*masse_so4
239			!JLD . /n_avogadro*1.e+12
240			so4_out(i,k) = sulfate(i,k)
241			IF (so4_out(i,k).LT.0) &
242			stop 'WAS SOLL DER SCHEISS ? '
243			ENDDO
244			ENDDO
245			ELSE ! if no new day, use old data:
246			DO k=1,klev
247			DO i=1,klon
248			sulfate(i,k) = so4_out(i,k)
249			IF (so4_out(i,k).LT.0) &
250			stop 'WAS SOLL DER SCHEISS ? '
251			ENDDO
252			ENDDO
253			ENDIF ! Did I have to do anything (was it a new day?)
254
255			END SUBROUTINE readsulfate
256
257			end module readsulfate_m