Sources/phylmd/readsulfate_preind.f

module readsulfate_preind_m

  IMPLICIT none

contains

  SUBROUTINE readsulfate_preind(dayvrai, time, first, pi_sulfate)

    ! Read in /calculate pre-industrial values of sulfate

    use dimens_m
    use dimphy
    use dynetat0_m, only: annee_ref
    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:
    ! ------
    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 pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data,
    !  from file)
    !
    ! Local Variables:
    ! ----------------
    INTEGER i, ig, k, it
    INTEGER j, iday, iyr

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


       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	130	SUBROUTINE readsulfate_preind(dayvrai, time, 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	guez	129	use dynetat0_m, only: annee_ref
14	guez	108	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	guez	130	integer, intent(in):: dayvrai
36			! current day number, based at value 1 on January 1st of annee_ref
37
38			REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
39
40	guez	108	LOGICAL, intent(in):: first ! First timestep
41			! (and therefore initialization necessary)
42			!
43			! Output:
44			! -------
45			real pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
46			! from file)
47			!
48			! Local Variables:
49			! ----------------
50			INTEGER i, ig, k, it
51	guez	188	INTEGER j, iday, iyr
52	guez	68
53	guez	108	INTEGER im, day1, day2, im2
54			double precision pi_so4_1(iim, jjm+1, klev, 12)
55	guez	68
56	guez	108	double precision pi_so4(klon, klev, 12) ! SO4 in right dimension
57			SAVE pi_so4
58			double precision pi_so4_out(klon, klev)
59			SAVE pi_so4_out
60	guez	68
61	guez	108	CHARACTER(len=4) cyear
62			LOGICAL lnewday
63	guez	68
64
65
66	guez	130	iday = dayvrai
67	guez	68
68	guez	108	! Get the year of the run
69			iyr = iday/360
70	guez	68
71	guez	108	! Get the day of the actual year:
72			iday = iday-iyr*360
73	guez	68
74	guez	108	! 0.02 is about 0.5/24, namly less than half an hour
75	guez	130	lnewday = time < 0.02
76	guez	68
77	guez	108	! ---------------------------------------------
78			! All has to be done only, if a new day begins!
79			! ---------------------------------------------
80	guez	68
81	guez	108	IF (lnewday.OR.first) THEN
82			im = iday/30 +1 ! the actual month
83	guez	68
84	guez	108	! annee_ref is the initial year (defined in temps.h)
85			iyr = iyr + annee_ref
86	guez	68
87
88	guez	108	IF (first) THEN
89			cyear='.nat'
90			CALL getso4fromfile(cyear,pi_so4_1)
91	guez	68
92	guez	108	! Transform the horizontal 2D-field into the physics-field
93			! (Also the levels and the latitudes have to be inversed)
94	guez	68
95	guez	108	! Initialize field
96			DO it=1,12
97			DO k=1,klev
98			DO i=1,klon
99			pi_so4(i,k,it)=0.
100			ENDDO
101			ENDDO
102			ENDDO
103	guez	68
104	guez	108	write (,) 'preind: finished reading', FLOAT(iim)
105			DO it=1,12
106			DO k=1, klev
107			! a) at the poles, use the zonal mean:
108			DO i=1,iim
109			! North pole
110			pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
111			! South pole
112			pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
113			ENDDO
114			pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
115			pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
116	guez	68
117	guez	108	! b) the values between the poles:
118			ig=1
119			DO j=2,jjm
120			DO i=1,iim
121			ig=ig+1
122			if (ig.gt.klon) write (,) 'shit'
123			pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)
124			ENDDO
125			ENDDO
126			IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
127			ENDDO ! Loop over k (vertical)
128			ENDDO ! Loop over it (months)
129	guez	68
130	guez	108	ENDIF ! Had to read new data?
131	guez	68
132
133	guez	108	! Interpolate to actual day:
134			IF (iday.LT.im*30-15) THEN
135			! in the first half of the month use month before and actual month
136			im2=im-1
137			day1 = im2*30+15
138			day2 = im2*30-15
139			IF (im2.LE.0) THEN
140			! the month is january, thus the month before december
141			im2=12
142			ENDIF
143			DO k=1,klev
144			DO i=1,klon
145			pi_sulfate(i,k) = pi_so4(i,k,im2) &
146			- FLOAT(iday-day2)/FLOAT(day1-day2) &
147			* (pi_so4(i,k,im2) - pi_so4(i,k,im))
148			IF (pi_sulfate(i,k).LT.0.) THEN
149			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
150			IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
151			write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
152			pi_so4(i,k,im2) - pi_so4(i,k,im)
153			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
154			stop 'pi_sulfate'
155			endif
156			ENDDO
157			ENDDO
158			ELSE
159			! the second half of the month
160			im2=im+1
161			day1 = im*30+15
162			IF (im2.GT.12) THEN
163			! the month is december, the following thus january
164			im2=1
165			ENDIF
166			day2 = im*30-15
167	guez	68
168	guez	108	DO k=1,klev
169			DO i=1,klon
170			pi_sulfate(i,k) = pi_so4(i,k,im2) &
171			- FLOAT(iday-day2)/FLOAT(day1-day2) &
172			* (pi_so4(i,k,im2) - pi_so4(i,k,im))
173			IF (pi_sulfate(i,k).LT.0.) THEN
174			IF (iday-day2.LT.0.) write(,) 'iday-day2',iday-day2
175			IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) &
176			write(,) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', &
177			pi_so4(i,k,im2) - pi_so4(i,k,im)
178			IF (day1-day2.LT.0.) write(,) 'day1-day2',day1-day2
179			stop 'pi_sulfate'
180			endif
181			ENDDO
182			ENDDO
183			ENDIF
184
185
186			!JLD ! The sulfate concentration [molec cm-3] is read in.
187			!JLD ! Convert it into mass [ug SO4/m3]
188			!JLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
189			DO k=1,klev
190			DO i=1,klon
191			!JLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
192			!JLD . /n_avogadro*1.e+12
193			pi_so4_out(i,k) = pi_sulfate(i,k)
194			ENDDO
195			ENDDO
196
197			ELSE ! If no new day, use old data:
198			DO k=1,klev
199			DO i=1,klon
200			pi_sulfate(i,k) = pi_so4_out(i,k)
201			ENDDO
202			ENDDO
203			ENDIF ! Was this the beginning of a new day?
204
205			END SUBROUTINE readsulfate_preind
206
207			end module readsulfate_preind_m