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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (show annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
File size: 17334 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21