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