/[lmdze]/trunk/Sources/phylmd/readsulfate.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/readsulfate.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 68 - (hide annotations)
Wed Nov 14 16:59:30 2012 UTC (11 years, 6 months ago) by guez
Original Path: trunk/libf/phylmd/readsulfate.f90
File size: 8073 byte(s) 
Split "flincom.f90" into "flinclo.f90", "flinfindcood.f90",
"flininfo.f90" and "flinopen_nozoom.f90", in directory
"IOIPSL/Flincom".

Renamed "etat0_lim" to "ce0l", as in LMDZ.

Split "readsulfate.f" into "readsulfate.f90", "readsulfate_preind.f90"
and "getso4fromfile.f90".

In etat0, renamed variable q3d to q, as in "dynredem1". Replaced calls
to Flicom procedures by calls to NetCDF95.

In leapfrog, added call to writehist.

Extracted ASCII art from "grid_noro" into a file
"grid_noro.txt". Transformed explicit-shape local arrays into
automatic arrays, so that test on values of iim and jjm is no longer
needed. Test on weight:
          IF (weight(ii, jj) /= 0.) THEN
is useless. There is already a test before:
    if (any(weight == 0.)) stop "zero weight in grid_noro"

In "aeropt", replaced duplicated lines with different values of inu by
a loop on inu.

Removed arguments of "conf_phys". Corresponding variables are now
defined in "physiq", in a namelist. In "conf_phys", read a namelist
instead of using getin.
1 guez 68 SUBROUTINE readsulfate(r_day, first, sulfate)
2 guez 3
3 guez 68 ! From LMDZ4/libf/phylmd/readsulfate.F,v 1.2 2005/05/19 08:27:15 fairhead
4 guez 3
5 guez 68 use dimens_m
6     use dimphy
7     use temps
8     use SUPHEC_M
9     use chem
10     IMPLICIT none
11 guez 3
12 guez 68 ! Content:
13     ! --------
14     ! This routine reads in monthly mean values of sulfate aerosols and
15     ! returns a linearly interpolated daily-mean field.
16     !
17     !
18     ! Author:
19     ! -------
20     ! Johannes Quaas (quaas@lmd.jussieu.fr)
21     ! 26/04/01
22     !
23     ! Modifications:
24     ! --------------
25     ! 21/06/01: Make integrations of more than one year possible ;-)
26     ! ATTENTION!! runs are supposed to start with Jan, 1. 1930
27     ! (rday=1)
28     !
29     ! 27/06/01: Correction: The model always has 360 days per year!
30     ! 27/06/01: SO4 concentration rather than mixing ratio
31     ! 27/06/01: 10yr-mean-values to interpolate
32     ! 20/08/01: Correct the error through integer-values in interpolations
33     ! 21/08/01: Introduce flag to read in just one decade
34     !
35     !
36     ! Input:
37     ! ------
38     double precision, intent(in):: r_day ! Day of integration
39     LOGICAL, intent(in):: first ! First timestep
40     ! (and therefore initialization necessary)
41     !
42     ! Output:
43     ! -------
44     double precision sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
45     ! from file) [ug SO4/m3]
46     !
47     ! Local Variables:
48     ! ----------------
49     INTEGER i, ig, k, it
50     INTEGER j, iday, ny, iyr, iyr1, iyr2
51     parameter (ny=jjm+1)
52 guez 3
53 guez 68 INTEGER ismaller
54     !JLD INTEGER idec1, idec2 ! The two decadal data read ini
55     CHARACTER*4 cyear
56 guez 3
57 guez 68 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 68 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 68 LOGICAL lnewday
67     LOGICAL lonlyone
68     PARAMETER (lonlyone=.FALSE.)
69 guez 3
70 guez 68 iday = INT(r_day)
71 guez 3
72 guez 68 ! Get the year of the run
73     iyr = iday/360
74 guez 3
75 guez 68 ! Get the day of the actual year:
76     iday = iday-iyr*360
77 guez 3
78 guez 68 ! 0.02 is about 0.5/24, namly less than half an hour
79     lnewday = (r_day-FLOAT(iday).LT.0.02)
80 guez 3
81 guez 68 ! ---------------------------------------------
82     ! All has to be done only, if a new day begins!
83     ! ---------------------------------------------
84 guez 3
85 guez 68 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 68 ! 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    
102 guez 68 IF (iyr .lt. 1850) THEN
103     cyear='.nat'
104     WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
105     CALL getso4fromfile(cyear, so4_1)
106     ELSE IF (iyr .ge. 2100) THEN
107     cyear='2100'
108     WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
109     CALL getso4fromfile(cyear, so4_1)
110     ELSE
111 guez 3
112 guez 68 ! Read in data:
113     ! a) from actual 10-yr-period
114 guez 3
115 guez 68 IF (iyr.LT.1900) THEN
116     iyr1 = 1850
117     iyr2 = 1900
118     ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
119     iyr1 = 1900
120     iyr2 = 1920
121     ELSE
122     iyr1 = INT(iyr/10)*10
123     iyr2 = INT(1+iyr/10)*10
124     ENDIF
125     WRITE(cyear,'(I4)') iyr1
126     WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
127     CALL getso4fromfile(cyear, so4_1)
128 guez 3
129    
130 guez 68 ! If to read two decades:
131     IF (.NOT.lonlyone) THEN
132 guez 3
133 guez 68 ! b) from the next following one
134     WRITE(cyear,'(I4)') iyr2
135     WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
136     CALL getso4fromfile(cyear, so4_2)
137 guez 3
138 guez 68 ENDIF
139 guez 3
140 guez 68 ! Interpolate linarily to the actual year:
141     DO it=1,12
142     DO k=1,klev
143     DO j=1,jjm
144     DO i=1,iim
145     so4_1(i,j,k,it)=so4_1(i,j,k,it) &
146     - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &
147     * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
148     ENDDO
149     ENDDO
150     ENDDO
151     ENDDO
152 guez 3
153 guez 68 ENDIF !lonlyone
154 guez 3
155 guez 68 ! Transform the horizontal 2D-field into the physics-field
156     ! (Also the levels and the latitudes have to be inversed)
157 guez 3
158 guez 68 DO it=1,12
159     DO k=1, klev
160     ! a) at the poles, use the zonal mean:
161     DO i=1,iim
162     ! North pole
163     so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
164     ! South pole
165     so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
166     ENDDO
167     so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
168     so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
169 guez 3
170 guez 68 ! b) the values between the poles:
171     ig=1
172     DO j=2,jjm
173     DO i=1,iim
174     ig=ig+1
175     if (ig.gt.klon) write (*,*) 'shit'
176     so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
177     ENDDO
178     ENDDO
179     IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
180     ENDDO ! Loop over k (vertical)
181     ENDDO ! Loop over it (months)
182 guez 3
183    
184 guez 68 ENDIF ! Had to read new data?
185 guez 3
186    
187 guez 68 ! Interpolate to actual day:
188     IF (iday.LT.im*30-15) THEN
189     ! in the first half of the month use month before and actual month
190     im2=im-1
191     day2 = im2*30-15
192     day1 = im2*30+15
193     IF (im2.LE.0) THEN
194     ! the month is january, thus the month before december
195     im2=12
196     ENDIF
197     DO k=1,klev
198     DO i=1,klon
199     sulfate(i,k) = so4(i,k,im2) &
200     - FLOAT(iday-day2)/FLOAT(day1-day2) &
201     * (so4(i,k,im2) - so4(i,k,im))
202     IF (sulfate(i,k).LT.0.) THEN
203     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
204     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
205     write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
206     so4(i,k,im2) - so4(i,k,im)
207     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
208     stop 'sulfate'
209     endif
210     ENDDO
211     ENDDO
212     ELSE
213     ! the second half of the month
214     im2=im+1
215     IF (im2.GT.12) THEN
216     ! the month is december, the following thus january
217     im2=1
218     ENDIF
219     day2 = im*30-15
220     day1 = im*30+15
221     DO k=1,klev
222     DO i=1,klon
223     sulfate(i,k) = so4(i,k,im2) &
224     - FLOAT(iday-day2)/FLOAT(day1-day2) &
225     * (so4(i,k,im2) - so4(i,k,im))
226     IF (sulfate(i,k).LT.0.) THEN
227     IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
228     IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
229     write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
230     so4(i,k,im2) - so4(i,k,im)
231     IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
232     stop 'sulfate'
233     endif
234     ENDDO
235     ENDDO
236     ENDIF
237 guez 3
238    
239 guez 68 !JLD ! The sulfate concentration [molec cm-3] is read in.
240     !JLD ! Convert it into mass [ug SO4/m3]
241     !JLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
242     ! The sulfate mass [ug SO4/m3] is read in.
243     DO k=1,klev
244     DO i=1,klon
245     !JLD sulfate(i,k) = sulfate(i,k)*masse_so4
246     !JLD . /n_avogadro*1.e+12
247     so4_out(i,k) = sulfate(i,k)
248     IF (so4_out(i,k).LT.0) &
249     stop 'WAS SOLL DER SCHEISS ? '
250     ENDDO
251     ENDDO
252     ELSE ! if no new day, use old data:
253     DO k=1,klev
254     DO i=1,klon
255     sulfate(i,k) = so4_out(i,k)
256     IF (so4_out(i,k).LT.0) &
257     stop 'WAS SOLL DER SCHEISS ? '
258     ENDDO
259     ENDDO
260     ENDIF ! Did I have to do anything (was it a new day?)
261 guez 3
262 guez 68 END SUBROUTINE readsulfate
  ViewVC Help
Powered by ViewVC 1.1.21