Context Navigation

readsulfate.F @ 1042

Last change on this file since 1042 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.8 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	#include "netcdf.inc"
540	#include "dimensions.h"
541	#include "dimphy.h"
542	CHARACTER*15 fname
543	CHARACTER*4 cyr
544
545	CHARACTER*6 cvar
546	INTEGER START(3), COUNT(3)
547	INTEGER STATUS, NCID, VARID
548	INTEGER imth, i, j, k, ny
549	PARAMETER (ny=jjm+1)
550
551
552	REAL so4mth(iim, ny, klev)
553	REAL so4(iim, ny, klev, 12)
554
555
556	fname = 'so4.run'//cyr//'.cdf'
557
558	write (,) 'reading ', fname
559	STATUS = NF_OPEN (fname, NF_NOWRITE, NCID)
560	IF (STATUS .NE. NF_NOERR) write (,) 'err in open ',status
561
562	DO imth=1, 12
563	IF (imth.eq.1) THEN
564	cvar='SO4JAN'
565	ELSEIF (imth.eq.2) THEN
566	cvar='SO4FEB'
567	ELSEIF (imth.eq.3) THEN
568	cvar='SO4MAR'
569	ELSEIF (imth.eq.4) THEN
570	cvar='SO4APR'
571	ELSEIF (imth.eq.5) THEN
572	cvar='SO4MAY'
573	ELSEIF (imth.eq.6) THEN
574	cvar='SO4JUN'
575	ELSEIF (imth.eq.7) THEN
576	cvar='SO4JUL'
577	ELSEIF (imth.eq.8) THEN
578	cvar='SO4AUG'
579	ELSEIF (imth.eq.9) THEN
580	cvar='SO4SEP'
581	ELSEIF (imth.eq.10) THEN
582	cvar='SO4OCT'
583	ELSEIF (imth.eq.11) THEN
584	cvar='SO4NOV'
585	ELSEIF (imth.eq.12) THEN
586	cvar='SO4DEC'
587	ENDIF
588	start(1)=1
589	start(2)=1
590	start(3)=1
591	count(1)=iim
592	count(2)=ny
593	count(3)=klev
594	c write(,) 'here i am'
595	STATUS = NF_INQ_VARID (NCID, cvar, VARID)
596	write (,) ncid,imth,cvar, varid
597
598	IF (STATUS .NE. NF_NOERR) write (,) 'err in read ',status
599
600	#ifdef NC_DOUBLE
601	status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth)
602	#else
603	status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth)
604	#endif
605	IF (STATUS .NE. NF_NOERR) write (,) 'err in read data',status
606
607	DO k=1,klev
608	DO j=1,jjm+1
609	DO i=1,iim
610	IF (so4mth(i,j,k).LT.0.) then
611	write(,) 'this is shit'
612	write(,) 'so4(',i,j,k,') =',so4mth(i,j,k)
613	endif
614	so4(i,j,k,imth)=so4mth(i,j,k)
615	ENDDO
616	ENDDO
617	ENDDO
618	ENDDO
619
620	STATUS = NF_CLOSE(NCID)
621	IF (STATUS .NE. NF_NOERR) write (,) 'err in closing file',status
622
623
624	END ! subroutine getso4fromfile
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640

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

Download in other formats:

Original Format