Context Navigation

readsulfate.F @ 857

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

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

Download in other formats:

Original Format