Sources/phylmd/readsulfate_preind.f

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
  IMPLICIT none

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