Context Navigation

readsulfate.F @ 1046

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

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

Download in other formats:

Original Format