Context Navigation

readsulfate.F @ 1261

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

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

Download in other formats:

Original Format