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

Contents of /trunk/libf/phylmd/readsulfate.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (show annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years ago) by guez
File size: 17236 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/readsulfate.F,v 1.2 2005/05/19 08:27:15 fairhead Exp $
3 !
4 SUBROUTINE readsulfate(r_day, first, sulfate)
5
6 use dimens_m
7 use dimphy
8 use temps
9 use SUPHEC_M
10 use chem
11 IMPLICIT none
12
13 c Content:
14 c --------
15 c This routine reads in monthly mean values of sulfate aerosols and
16 c returns a linearly interpolated dayly-mean field.
17 c
18 c
19 c Author:
20 c -------
21 c Johannes Quaas (quaas@lmd.jussieu.fr)
22 c 26/04/01
23 c
24 c Modifications:
25 c --------------
26 c 21/06/01: Make integrations of more than one year possible ;-)
27 c ATTENTION!! runs are supposed to start with Jan, 1. 1930
28 c (rday=1)
29 c
30 c 27/06/01: Correction: The model always has 360 days per year!
31 c 27/06/01: SO4 concentration rather than mixing ratio
32 c 27/06/01: 10yr-mean-values to interpolate
33 c 20/08/01: Correct the error through integer-values in interpolations
34 c 21/08/01: Introduce flag to read in just one decade
35 c
36 c
37 c Input:
38 c ------
39 REAL*8, intent(in):: r_day ! Day of integration
40 LOGICAL, intent(in):: first ! First timestep
41 ! (and therefore initialization necessary)
42 c
43 c Output:
44 c -------
45 REAL*8 sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
46 ! from file) [ug SO4/m3]
47 c
48 c Local Variables:
49 c ----------------
50 INTEGER i, ig, k, it
51 INTEGER j, iday, ny, iyr, iyr1, iyr2
52 parameter (ny=jjm+1)
53
54 INTEGER ismaller
55 CJLD INTEGER idec1, idec2 ! The two decadal data read ini
56 CHARACTER*4 cyear
57
58 INTEGER im, day1, day2, im2
59 REAL*8 so4_1(iim, jjm+1, klev, 12)
60 REAL*8 so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions
61
62 REAL*8 so4(klon, klev, 12) ! SO4 in right dimension
63 SAVE so4
64 REAL*8 so4_out(klon, klev)
65 SAVE so4_out
66
67 LOGICAL lnewday
68 LOGICAL lonlyone
69 PARAMETER (lonlyone=.FALSE.)
70
71 iday = INT(r_day)
72
73 ! Get the year of the run
74 iyr = iday/360
75
76 ! Get the day of the actual year:
77 iday = iday-iyr*360
78
79 ! 0.02 is about 0.5/24, namly less than half an hour
80 lnewday = (r_day-FLOAT(iday).LT.0.02)
81
82 ! ---------------------------------------------
83 ! All has to be done only, if a new day begins!
84 ! ---------------------------------------------
85
86 IF (lnewday.OR.first) THEN
87
88 im = iday/30 +1 ! the actual month
89 ! annee_ref is the initial year (defined in temps.h)
90 iyr = iyr + annee_ref
91
92 ! Do I have to read new data? (Is this the first day of a year?)
93 IF (first.OR.iday.EQ.1.) THEN
94 ! Initialize values
95 DO it=1,12
96 DO k=1,klev
97 DO i=1,klon
98 so4(i,k,it)=0.
99 ENDDO
100 ENDDO
101 ENDDO
102
103
104 IF (iyr .lt. 1850) THEN
105 cyear='.nat'
106 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
107 CALL getso4fromfile(cyear, so4_1)
108 ELSE IF (iyr .ge. 2100) THEN
109 cyear='2100'
110 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
111 CALL getso4fromfile(cyear, so4_1)
112 ELSE
113
114 ! Read in data:
115 ! a) from actual 10-yr-period
116
117 IF (iyr.LT.1900) THEN
118 iyr1 = 1850
119 iyr2 = 1900
120 ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
121 iyr1 = 1900
122 iyr2 = 1920
123 ELSE
124 iyr1 = INT(iyr/10)*10
125 iyr2 = INT(1+iyr/10)*10
126 ENDIF
127 WRITE(cyear,'(I4)') iyr1
128 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
129 CALL getso4fromfile(cyear, so4_1)
130
131
132 ! If to read two decades:
133 IF (.NOT.lonlyone) THEN
134
135 ! b) from the next following one
136 WRITE(cyear,'(I4)') iyr2
137 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear
138 CALL getso4fromfile(cyear, so4_2)
139
140 ENDIF
141
142 ! Interpolate linarily to the actual year:
143 DO it=1,12
144 DO k=1,klev
145 DO j=1,jjm
146 DO i=1,iim
147 so4_1(i,j,k,it)=so4_1(i,j,k,it)
148 . - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
149 . * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
150 ENDDO
151 ENDDO
152 ENDDO
153 ENDDO
154
155 ENDIF !lonlyone
156
157 ! Transform the horizontal 2D-field into the physics-field
158 ! (Also the levels and the latitudes have to be inversed)
159
160 DO it=1,12
161 DO k=1, klev
162 ! a) at the poles, use the zonal mean:
163 DO i=1,iim
164 ! North pole
165 so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
166 ! South pole
167 so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
168 ENDDO
169 so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
170 so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
171
172 ! b) the values between the poles:
173 ig=1
174 DO j=2,jjm
175 DO i=1,iim
176 ig=ig+1
177 if (ig.gt.klon) write (*,*) 'shit'
178 so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
179 ENDDO
180 ENDDO
181 IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
182 ENDDO ! Loop over k (vertical)
183 ENDDO ! Loop over it (months)
184
185
186 ENDIF ! Had to read new data?
187
188
189 ! Interpolate to actual day:
190 IF (iday.LT.im*30-15) THEN
191 ! in the first half of the month use month before and actual month
192 im2=im-1
193 day2 = im2*30-15
194 day1 = im2*30+15
195 IF (im2.LE.0) THEN
196 ! the month is january, thus the month before december
197 im2=12
198 ENDIF
199 DO k=1,klev
200 DO i=1,klon
201 sulfate(i,k) = so4(i,k,im2)
202 . - FLOAT(iday-day2)/FLOAT(day1-day2)
203 . * (so4(i,k,im2) - so4(i,k,im))
204 IF (sulfate(i,k).LT.0.) THEN
205 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
206 IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
207 . write(*,*) 'so4(i,k,im2) - so4(i,k,im)',
208 . so4(i,k,im2) - so4(i,k,im)
209 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
210 stop 'sulfate'
211 endif
212 ENDDO
213 ENDDO
214 ELSE
215 ! the second half of the month
216 im2=im+1
217 IF (im2.GT.12) THEN
218 ! the month is december, the following thus january
219 im2=1
220 ENDIF
221 day2 = im*30-15
222 day1 = im*30+15
223 DO k=1,klev
224 DO i=1,klon
225 sulfate(i,k) = so4(i,k,im2)
226 . - FLOAT(iday-day2)/FLOAT(day1-day2)
227 . * (so4(i,k,im2) - so4(i,k,im))
228 IF (sulfate(i,k).LT.0.) THEN
229 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
230 IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
231 . write(*,*) 'so4(i,k,im2) - so4(i,k,im)',
232 . so4(i,k,im2) - so4(i,k,im)
233 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
234 stop 'sulfate'
235 endif
236 ENDDO
237 ENDDO
238 ENDIF
239
240
241 CJLD ! The sulfate concentration [molec cm-3] is read in.
242 CJLD ! Convert it into mass [ug SO4/m3]
243 CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
244 ! The sulfate mass [ug SO4/m3] is read in.
245 DO k=1,klev
246 DO i=1,klon
247 CJLD sulfate(i,k) = sulfate(i,k)*masse_so4
248 CJLD . /n_avogadro*1.e+12
249 so4_out(i,k) = sulfate(i,k)
250 IF (so4_out(i,k).LT.0)
251 . stop 'WAS SOLL DER SCHEISS ? '
252 ENDDO
253 ENDDO
254 ELSE ! if no new day, use old data:
255 DO k=1,klev
256 DO i=1,klon
257 sulfate(i,k) = so4_out(i,k)
258 IF (so4_out(i,k).LT.0)
259 . stop 'WAS SOLL DER SCHEISS ? '
260 ENDDO
261 ENDDO
262
263
264 ENDIF ! Did I have to do anything (was it a new day?)
265
266 RETURN
267 END
268
269
270
271
272
273 c-----------------------------------------------------------------------------
274 c Read in /calculate pre-industrial values of sulfate
275 c-----------------------------------------------------------------------------
276
277 SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate)
278
279 use dimens_m
280 use dimphy
281 use temps
282 use SUPHEC_M
283 use chem
284 IMPLICIT none
285
286 c Content:
287 c --------
288 c This routine reads in monthly mean values of sulfate aerosols and
289 c returns a linearly interpolated dayly-mean field.
290 c
291 c It does so for the preindustriel values of the sulfate, to a large part
292 c analogous to the routine readsulfate above.
293 c
294 c Only Pb: Variables must be saved and don t have to be overwritten!
295 c
296 c Author:
297 c -------
298 c Johannes Quaas (quaas@lmd.jussieu.fr)
299 c 26/06/01
300 c
301 c Modifications:
302 c --------------
303 c see above
304 c
305 c
306 c Input:
307 c ------
308 REAL*8, intent(in):: r_day ! Day of integration
309 LOGICAL, intent(in):: first ! First timestep
310 ! (and therefore initialization necessary)
311 c
312 c Output:
313 c -------
314 REAL*8 pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
315 ! from file)
316 c
317 c Local Variables:
318 c ----------------
319 INTEGER i, ig, k, it
320 INTEGER j, iday, ny, iyr
321 parameter (ny=jjm+1)
322
323 INTEGER im, day1, day2, im2, ismaller
324 REAL*8 pi_so4_1(iim, jjm+1, klev, 12)
325
326 REAL*8 pi_so4(klon, klev, 12) ! SO4 in right dimension
327 SAVE pi_so4
328 REAL*8 pi_so4_out(klon, klev)
329 SAVE pi_so4_out
330
331 CHARACTER*4 cyear
332 LOGICAL lnewday
333
334
335
336 iday = INT(r_day)
337
338 ! Get the year of the run
339 iyr = iday/360
340
341 ! Get the day of the actual year:
342 iday = iday-iyr*360
343
344 ! 0.02 is about 0.5/24, namly less than half an hour
345 lnewday = (r_day-FLOAT(iday).LT.0.02)
346
347 ! ---------------------------------------------
348 ! All has to be done only, if a new day begins!
349 ! ---------------------------------------------
350
351 IF (lnewday.OR.first) THEN
352
353 im = iday/30 +1 ! the actual month
354
355 ! annee_ref is the initial year (defined in temps.h)
356 iyr = iyr + annee_ref
357
358
359 IF (first) THEN
360 cyear='.nat'
361 CALL getso4fromfile(cyear,pi_so4_1)
362
363 ! Transform the horizontal 2D-field into the physics-field
364 ! (Also the levels and the latitudes have to be inversed)
365
366 ! Initialize field
367 DO it=1,12
368 DO k=1,klev
369 DO i=1,klon
370 pi_so4(i,k,it)=0.
371 ENDDO
372 ENDDO
373 ENDDO
374
375 write (*,*) 'preind: finished reading', FLOAT(iim)
376 DO it=1,12
377 DO k=1, klev
378 ! a) at the poles, use the zonal mean:
379 DO i=1,iim
380 ! North pole
381 pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
382 ! South pole
383 pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
384 ENDDO
385 pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
386 pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
387
388 ! b) the values between the poles:
389 ig=1
390 DO j=2,jjm
391 DO i=1,iim
392 ig=ig+1
393 if (ig.gt.klon) write (*,*) 'shit'
394 pi_so4(ig,k,it) = pi_so4_1(i,jjm+1-j,klev+1-k,it)
395 ENDDO
396 ENDDO
397 IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
398 ENDDO ! Loop over k (vertical)
399 ENDDO ! Loop over it (months)
400
401 ENDIF ! Had to read new data?
402
403
404 ! Interpolate to actual day:
405 IF (iday.LT.im*30-15) THEN
406 ! in the first half of the month use month before and actual month
407 im2=im-1
408 day1 = im2*30+15
409 day2 = im2*30-15
410 IF (im2.LE.0) THEN
411 ! the month is january, thus the month before december
412 im2=12
413 ENDIF
414 DO k=1,klev
415 DO i=1,klon
416 pi_sulfate(i,k) = pi_so4(i,k,im2)
417 . - FLOAT(iday-day2)/FLOAT(day1-day2)
418 . * (pi_so4(i,k,im2) - pi_so4(i,k,im))
419 IF (pi_sulfate(i,k).LT.0.) THEN
420 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
421 IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
422 . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
423 . pi_so4(i,k,im2) - pi_so4(i,k,im)
424 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
425 stop 'pi_sulfate'
426 endif
427 ENDDO
428 ENDDO
429 ELSE
430 ! the second half of the month
431 im2=im+1
432 day1 = im*30+15
433 IF (im2.GT.12) THEN
434 ! the month is december, the following thus january
435 im2=1
436 ENDIF
437 day2 = im*30-15
438
439 DO k=1,klev
440 DO i=1,klon
441 pi_sulfate(i,k) = pi_so4(i,k,im2)
442 . - FLOAT(iday-day2)/FLOAT(day1-day2)
443 . * (pi_so4(i,k,im2) - pi_so4(i,k,im))
444 IF (pi_sulfate(i,k).LT.0.) THEN
445 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
446 IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
447 . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
448 . pi_so4(i,k,im2) - pi_so4(i,k,im)
449 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
450 stop 'pi_sulfate'
451 endif
452 ENDDO
453 ENDDO
454 ENDIF
455
456
457 CJLD ! The sulfate concentration [molec cm-3] is read in.
458 CJLD ! Convert it into mass [ug SO4/m3]
459 CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
460 DO k=1,klev
461 DO i=1,klon
462 CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
463 CJLD . /n_avogadro*1.e+12
464 pi_so4_out(i,k) = pi_sulfate(i,k)
465 ENDDO
466 ENDDO
467
468 ELSE ! If no new day, use old data:
469 DO k=1,klev
470 DO i=1,klon
471 pi_sulfate(i,k) = pi_so4_out(i,k)
472 ENDDO
473 ENDDO
474
475
476 ENDIF ! Was this the beginning of a new day?
477 RETURN
478 END
479
480
481
482
483
484
485
486
487
488
489 c-----------------------------------------------------------------------------
490 c Routine for reading SO4 data from files
491 c-----------------------------------------------------------------------------
492
493
494 SUBROUTINE getso4fromfile (cyr, so4)
495 use dimens_m
496 use dimphy
497 use netcdf
498 CHARACTER*15 fname
499 CHARACTER*4 cyr
500
501 CHARACTER*6 cvar
502 INTEGER START(3), COUNT(3)
503 INTEGER STATUS, NCID, VARID
504 INTEGER imth, i, j, k, ny
505 PARAMETER (ny=jjm+1)
506
507
508 REAL*8 so4mth(iim, ny, klev)
509 c REAL*8 so4mth(klev, ny, iim)
510 REAL*8 so4(iim, ny, klev, 12)
511
512
513 fname = 'so4.run'//cyr//'.cdf'
514
515 write (*,*) 'reading ', fname
516 STATUS = NF_OPEN (fname, NF_NOWRITE, NCID)
517 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in open ',status
518
519 DO imth=1, 12
520 IF (imth.eq.1) THEN
521 cvar='SO4JAN'
522 ELSEIF (imth.eq.2) THEN
523 cvar='SO4FEB'
524 ELSEIF (imth.eq.3) THEN
525 cvar='SO4MAR'
526 ELSEIF (imth.eq.4) THEN
527 cvar='SO4APR'
528 ELSEIF (imth.eq.5) THEN
529 cvar='SO4MAY'
530 ELSEIF (imth.eq.6) THEN
531 cvar='SO4JUN'
532 ELSEIF (imth.eq.7) THEN
533 cvar='SO4JUL'
534 ELSEIF (imth.eq.8) THEN
535 cvar='SO4AUG'
536 ELSEIF (imth.eq.9) THEN
537 cvar='SO4SEP'
538 ELSEIF (imth.eq.10) THEN
539 cvar='SO4OCT'
540 ELSEIF (imth.eq.11) THEN
541 cvar='SO4NOV'
542 ELSEIF (imth.eq.12) THEN
543 cvar='SO4DEC'
544 ENDIF
545 start(1)=1
546 start(2)=1
547 start(3)=1
548 count(1)=iim
549 count(2)=ny
550 count(3)=klev
551 c write(*,*) 'here i am'
552 STATUS = NF_INQ_VARID (NCID, cvar, VARID)
553 write (*,*) ncid,imth,cvar, varid
554 c STATUS = NF_INQ_VARID (NCID, VARMONTHS(i), VARID(i))
555 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read ',status
556 STATUS = NF_GET_VARA_DOUBLE
557 . (NCID, VARID, START,COUNT, so4mth)
558 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read data',status
559
560 DO k=1,klev
561 DO j=1,jjm+1
562 DO i=1,iim
563 IF (so4mth(i,j,k).LT.0.) then
564 write(*,*) 'this is shit'
565 write(*,*) 'so4(',i,j,k,') =',so4mth(i,j,k)
566 endif
567 so4(i,j,k,imth)=so4mth(i,j,k)
568 c so4(i,j,k,imth)=so4mth(k,j,i)
569 ENDDO
570 ENDDO
571 ENDDO
572 ENDDO
573
574 STATUS = NF_CLOSE(NCID)
575 END ! subroutine getso4fromfile
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591

  ViewVC Help
Powered by ViewVC 1.1.21