Context Navigation

readsulfate.F @ 1206

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

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

Download in other formats:

Original Format