Sources/phylmd/readsulfate.f

module readsulfate_m

  IMPLICIT none

contains

  SUBROUTINE readsulfate(dayvrai, time, 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 dynetat0_m, ONLY: annee_ref
    use getso4fromfile_m, only: getso4fromfile

    ! Input:

    integer, intent(in):: dayvrai
    ! current day number, based at value 1 on January 1st of annee_ref

    REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour

    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, iyr, iyr1, iyr2

    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 = dayvrai

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