Context Navigation

readsulfate.F @ 776

Last change on this file since 776 was 776, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format