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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 51 - (show annotations)
Tue Sep 20 09:14:34 2011 UTC (12 years, 7 months ago) by guez
File size: 17318 byte(s)
Split "getincom.f90" into "getincom.f90" and "getincom2.f90". Split
"nuage.f" into "nuage.f90", "diagcld1.f90" and "diagcld2.f90". Created
module "chem" from included file "chem.h". Moved "YOEGWD.f90" to
directory "Orography".

In "physiq", for evaporation of water, "zlsdcp" was equal to
"zlvdc". Removed useless variables.

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 Include-files:
37 c --------------
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 use chem
286 IMPLICIT none
287
288 c Content:
289 c --------
290 c This routine reads in monthly mean values of sulfate aerosols and
291 c returns a linearly interpolated dayly-mean field.
292 c
293 c It does so for the preindustriel values of the sulfate, to a large part
294 c analogous to the routine readsulfate above.
295 c
296 c Only Pb: Variables must be saved and don t have to be overwritten!
297 c
298 c Author:
299 c -------
300 c Johannes Quaas (quaas@lmd.jussieu.fr)
301 c 26/06/01
302 c
303 c Modifications:
304 c --------------
305 c see above
306 c
307 c Include-files:
308 c --------------
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