Sources/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)

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