Context Navigation

readsulfate.F @ 2192

Last change on this file since 2192 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.4 KB

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

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

Download in other formats:

Original Format