Sources/phylmd/readsulfate_preind.f

module readsulfate_preind_m

  IMPLICIT none

contains

  SUBROUTINE readsulfate_preind(r_day, first, pi_sulfate)

    ! 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

    ! Content:
    ! --------
    ! This routine reads in monthly mean values of sulfate aerosols and
    ! returns a linearly interpolated daily-mean field.
    !
    ! It does so for the preindustriel values of the sulfate, to a large part
    ! analogous to the routine readsulfate.
    !
    ! Only Pb: Variables must be saved and don t have to be overwritten!
    !
    ! Author:
    ! -------
    ! Johannes Quaas (quaas@lmd.jussieu.fr)
    ! 26/06/01
    !
    ! Input:
    ! ------
    real, intent(in)::  r_day                   ! Day of integration
    LOGICAL, intent(in):: first                 ! First timestep
    ! (and therefore initialization necessary)
    !
    ! Output:
    ! -------
    real pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data,
    !  from file)
    !
    ! Local Variables:
    ! ----------------
    INTEGER i, ig, k, it
    INTEGER j, iday, ny, iyr
    parameter (ny=jjm+1)

    INTEGER im, day1, day2, im2
    double precision pi_so4_1(iim, jjm+1, klev, 12)

    double precision pi_so4(klon, klev, 12)  ! SO4 in right dimension
    SAVE pi_so4
    double precision pi_so4_out(klon, klev)
    SAVE pi_so4_out

    CHARACTER(len=4) cyear
    LOGICAL lnewday


    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


       IF (first) THEN
          cyear='.nat'
          CALL getso4fromfile(cyear,pi_so4_1)

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

          ! Initialize field
          DO it=1,12
             DO k=1,klev
                DO i=1,klon
                   pi_so4(i,k,it)=0.
                ENDDO
             ENDDO
          ENDDO

          write (*,*) 'preind: finished reading', FLOAT(iim)
          DO it=1,12
             DO k=1, klev
                ! a) at the poles, use the zonal mean:
                DO i=1,iim
                   ! North pole
                   pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
                   ! South pole
                   pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
                ENDDO
                pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
                pi_so4(klon,k,it)=pi_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'
                      pi_so4(ig,k,it) = pi_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
          day1 = im2*30+15
          day2 = 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
                pi_sulfate(i,k) = pi_so4(i,k,im2)   &
                     - FLOAT(iday-day2)/FLOAT(day1-day2) &
                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))
                IF (pi_sulfate(i,k).LT.0.) THEN
                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
                        pi_so4(i,k,im2) - pi_so4(i,k,im)
                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
                   stop 'pi_sulfate'
                endif
             ENDDO
          ENDDO
       ELSE
          ! the second half of the month
          im2=im+1
          day1 = im*30+15
          IF (im2.GT.12) THEN
             ! the month is december, the following thus january
             im2=1
          ENDIF
          day2 = im*30-15

          DO k=1,klev
             DO i=1,klon
                pi_sulfate(i,k) = pi_so4(i,k,im2)   &
                     - FLOAT(iday-day2)/FLOAT(day1-day2) &
                     * (pi_so4(i,k,im2) - pi_so4(i,k,im))
                IF (pi_sulfate(i,k).LT.0.) THEN
                   IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
                   IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
                        write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
                        pi_so4(i,k,im2) - pi_so4(i,k,im)
                   IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
                   stop 'pi_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]
       DO k=1,klev
          DO i=1,klon
             !JLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
             !JLD     .           /n_avogadro*1.e+12
             pi_so4_out(i,k) = pi_sulfate(i,k)
          ENDDO
       ENDDO

    ELSE ! If no new day, use old data:
       DO k=1,klev
          DO i=1,klon
             pi_sulfate(i,k) = pi_so4_out(i,k)
          ENDDO
       ENDDO
    ENDIF ! Was this the beginning of a new day?

  END SUBROUTINE readsulfate_preind

end module readsulfate_preind_m
1	guez	108	module readsulfate_preind_m
2	guez	68
3			IMPLICIT none
4
5	guez	108	contains
6	guez	68
7	guez	108	SUBROUTINE readsulfate_preind(r_day, first, pi_sulfate)
8	guez	68
9	guez	108	! Read in /calculate pre-industrial values of sulfate
10	guez	68
11	guez	108	use dimens_m
12			use dimphy
13			use temps
14			use SUPHEC_M
15			use chem
16			use getso4fromfile_m, only: getso4fromfile
17	guez	68
18	guez	108	! Content:
19			! --------
20			! This routine reads in monthly mean values of sulfate aerosols and
21			! returns a linearly interpolated daily-mean field.
22			!
23			! It does so for the preindustriel values of the sulfate, to a large part
24			! analogous to the routine readsulfate.
25			!
26			! Only Pb: Variables must be saved and don t have to be overwritten!
27			!
28			! Author:
29			! -------
30			! Johannes Quaas (quaas@lmd.jussieu.fr)
31			! 26/06/01
32			!
33			! Input:
34			! ------
35			real, intent(in):: r_day ! Day of integration
36			LOGICAL, intent(in):: first ! First timestep
37			! (and therefore initialization necessary)
38			!
39			! Output:
40			! -------
41			real pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
42			! from file)
43			!
44			! Local Variables:
45			! ----------------
46			INTEGER i, ig, k, it
47			INTEGER j, iday, ny, iyr
48			parameter (ny=jjm+1)
49	guez	68
50	guez	108	INTEGER im, day1, day2, im2
51			double precision pi_so4_1(iim, jjm+1, klev, 12)
52	guez	68
53	guez	108	double precision pi_so4(klon, klev, 12) ! SO4 in right dimension
54			SAVE pi_so4
55			double precision pi_so4_out(klon, klev)
56			SAVE pi_so4_out
57	guez	68
58	guez	108	CHARACTER(len=4) cyear
59			LOGICAL lnewday
60	guez	68
61
62
63	guez	108	iday = INT(r_day)
64	guez	68
65	guez	108	! Get the year of the run
66			iyr = iday/360
67	guez	68
68	guez	108	! Get the day of the actual year:
69			iday = iday-iyr*360
70	guez	68
71	guez	108	! 0.02 is about 0.5/24, namly less than half an hour
72			lnewday = (r_day-FLOAT(iday).LT.0.02)
73	guez	68
74	guez	108	! ---------------------------------------------
75			! All has to be done only, if a new day begins!
76			! ---------------------------------------------
77	guez	68
78	guez	108	IF (lnewday.OR.first) THEN
79			im = iday/30 +1 ! the actual month
80	guez	68
81	guez	108	! annee_ref is the initial year (defined in temps.h)
82			iyr = iyr + annee_ref
83	guez	68
84
85	guez	108	IF (first) THEN
86			cyear='.nat'
87			CALL getso4fromfile(cyear,pi_so4_1)
88	guez	68
89	guez	108	! Transform the horizontal 2D-field into the physics-field
90			! (Also the levels and the latitudes have to be inversed)
91	guez	68
92	guez	108	! Initialize field
93			DO it=1,12
94			DO k=1,klev
95			DO i=1,klon
96			pi_so4(i,k,it)=0.
97			ENDDO
98			ENDDO
99			ENDDO
100	guez	68
101	guez	108	write (,) 'preind: finished reading', FLOAT(iim)
102			DO it=1,12
103			DO k=1, klev
104			! a) at the poles, use the zonal mean:
105			DO i=1,iim
106			! North pole
107			pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
108			! South pole
109			pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
110			ENDDO
111			pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
112			pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
113	guez	68
114	guez	108	! b) the values between the poles:
115			ig=1
116			DO j=2,jjm
117			DO i=1,iim
118			ig=ig+1
119			if (ig.gt.klon) write (,) 'shit'
120			pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)
121			ENDDO
122			ENDDO
123			IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
124			ENDDO ! Loop over k (vertical)
125			ENDDO ! Loop over it (months)
126	guez	68
127	guez	108	ENDIF ! Had to read new data?
128	guez	68
129
130	guez	108	! Interpolate to actual day:
131			IF (iday.LT.im*30-15) THEN
132			! in the first half of the month use month before and actual month
133			im2=im-1
134			day1 = im2*30+15
135			day2 = im2*30-15
136			IF (im2.LE.0) THEN
137			! the month is january, thus the month before december
138			im2=12
139			ENDIF
140			DO k=1,klev
141			DO i=1,klon
142			pi_sulfate(i,k) = pi_so4(i,k,im2) &
143			- FLOAT(iday-day2)/FLOAT(day1-day2) &
144			* (pi_so4(i,k,im2) - pi_so4(i,k,im))
145			IF (pi_sulfate(i,k).LT.0.) THEN
146			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
147			IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
148			write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
149			pi_so4(i,k,im2) - pi_so4(i,k,im)
150			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
151			stop 'pi_sulfate'
152			endif
153			ENDDO
154			ENDDO
155			ELSE
156			! the second half of the month
157			im2=im+1
158			day1 = im*30+15
159			IF (im2.GT.12) THEN
160			! the month is december, the following thus january
161			im2=1
162			ENDIF
163			day2 = im*30-15
164	guez	68
165	guez	108	DO k=1,klev
166			DO i=1,klon
167			pi_sulfate(i,k) = pi_so4(i,k,im2) &
168			- FLOAT(iday-day2)/FLOAT(day1-day2) &
169			* (pi_so4(i,k,im2) - pi_so4(i,k,im))
170			IF (pi_sulfate(i,k).LT.0.) THEN
171			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
172			IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
173			write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
174			pi_so4(i,k,im2) - pi_so4(i,k,im)
175			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
176			stop 'pi_sulfate'
177			endif
178			ENDDO
179			ENDDO
180			ENDIF
181
182
183			!JLD ! The sulfate concentration [molec cm-3] is read in.
184			!JLD ! Convert it into mass [ug SO4/m3]
185			!JLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
186			DO k=1,klev
187			DO i=1,klon
188			!JLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
189			!JLD . /n_avogadro*1.e+12
190			pi_so4_out(i,k) = pi_sulfate(i,k)
191			ENDDO
192			ENDDO
193
194			ELSE ! If no new day, use old data:
195			DO k=1,klev
196			DO i=1,klon
197			pi_sulfate(i,k) = pi_so4_out(i,k)
198			ENDDO
199			ENDDO
200			ENDIF ! Was this the beginning of a new day?
201
202			END SUBROUTINE readsulfate_preind
203
204			end module readsulfate_preind_m