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