Context Navigation

readsulfate.F @ 1077

Last change on this file since 1077 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

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
[766]	4	SUBROUTINE readsulfate (r_day, first, sulfate_p)
[776]	5	USE dimphy, ONLY : klev
	6	USE mod_grid_phy_lmdz, klon=>klon_glo
	7	USE mod_phys_lmdz_para
[524]	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"
[766]	38	cym#include "dimphy.h"
[524]	39	#include "temps.h"
	40	c
	41	c Input:
	42	c ------
[782]	43	REAL r_day ! Day of integration
[524]	44	LOGICAL first ! First timestep
	45	! (and therefore initialization necessary)
	46	c
	47	c Output:
	48	c -------
[782]	49	REAL sulfate_p(klon_omp,klev)
	50	REAL sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
[524]	51	! from file) [ug SO4/m3]
	52	c
	53	c Local Variables:
	54	c ----------------
	55	INTEGER i, ig, k, it
[640]	56	INTEGER j, iday, ny, iyr, iyr1, iyr2
[524]	57	parameter (ny=jjm+1)
	58
[640]	59	CJLD INTEGER idec1, idec2 ! The two decadal data read ini
[524]	60	CHARACTER*4 cyear
	61
	62	INTEGER im, day1, day2, im2
[782]	63	REAL so4_1(iim, jjm+1, klev, 12)
	64	REAL so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions
[524]	65
[782]	66	REAL, allocatable,save :: so4(:, :, :) ! SO4 in right dimension
	67	REAL, allocatable,save :: so4_out(:, :)
[766]	68	c$OMP THREADPRIVATE(so4,so4_out)
[524]	69
	70	LOGICAL lnewday
	71	LOGICAL lonlyone
	72	PARAMETER (lonlyone=.FALSE.)
[766]	73	logical,save :: first2=.true.
	74	c$OMP THREADPRIVATE(first2)
[524]	75
[766]	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
[776]	85	if (is_mpi_root) then
[766]	86
[524]	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
[640]	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:
[524]	131	! a) from actual 10-yr-period
	132
[640]	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
[524]	142	ENDIF
[640]	143	WRITE(cyear,'(I4)') iyr1
	144	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
[524]	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
[640]	152	WRITE(cyear,'(I4)') iyr2
	153	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
	154	CALL getso4fromfile(cyear, so4_2)
	155
[524]	156	ENDIF
[640]	157
[524]	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)
[640]	164	. - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
[524]	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'
[804]	194	so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	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
[640]	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.
[524]	261	DO k=1,klev
	262	DO i=1,klon
[640]	263	CJLD sulfate(i,k) = sulfate(i,k)*masse_so4
	264	CJLD . /n_avogadro*1.e+12
[524]	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?)
[766]	281
	282	endif ! phy_rank==0
[524]	283
[766]	284	c$OMP END MASTER
[782]	285	call Scatter(sulfate,sulfate_p)
[766]	286
[524]	287	RETURN
	288	END
	289
	290
	291
	292
	293
	294	c-----------------------------------------------------------------------------
	295	c Read in /calculate pre-industrial values of sulfate
	296	c-----------------------------------------------------------------------------
	297
[766]	298	SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p)
[776]	299	USE dimphy, ONLY : klev
	300	USE mod_grid_phy_lmdz, klon=>klon_glo
	301	USE mod_phys_lmdz_para
[524]	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 ------
[782]	332	REAL r_day ! Day of integration
[524]	333	LOGICAL first ! First timestep
	334	! (and therefore initialization necessary)
	335	c
	336	c Output:
	337	c -------
[782]	338	REAL pi_sulfate_p (klon_omp, klev)
[766]	339
[782]	340	REAL pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
[766]	341	! from fil
[524]	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
[782]	349	INTEGER im, day1, day2, im2
	350	REAL pi_so4_1(iim, jjm+1, klev, 12)
[766]	351
[782]	352	REAL, allocatable,save :: pi_so4(:, :, :) ! SO4 in right dimension
	353	REAL, allocatable,save :: pi_so4_out(:, :)
[766]	354	c$OMP THREADPRIVATE(pi_so4,pi_so4_out)
[524]	355
	356	CHARACTER*4 cyear
	357	LOGICAL lnewday
[766]	358	logical,save :: first2=.true.
	359	c$OMP THREADPRIVATE(first2)
[524]	360
[766]	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
[776]	370	if (is_mpi_root) then
[766]	371
[524]	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'
[804]	432	pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	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
[640]	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]
[524]	498	DO k=1,klev
	499	DO i=1,klon
[640]	500	CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
	501	CJLD . /n_avogadro*1.e+12
[524]	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?
[766]	515
[776]	516	endif ! is_mpi_root==0
[766]	517
	518	c$OMP END MASTER
[782]	519	call Scatter(pi_sulfate,pi_sulfate_p)
[766]	520
[524]	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)
[940]	539	use dimphy
[524]	540	#include "netcdf.inc"
	541	#include "dimensions.h"
[940]	542	cccc#include "dimphy.h"
[524]	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
[782]	553	REAL so4mth(iim, ny, klev)
	554	REAL so4(iim, ny, klev, 12)
[524]	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
[782]	598
[524]	599	IF (STATUS .NE. NF_NOERR) write (,) 'err in read ',status
[782]	600
	601	#ifdef NC_DOUBLE
[801]	602	status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth)
[782]	603	#else
[801]	604	status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth)
[782]	605	#endif
[524]	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)
[782]	622	IF (STATUS .NE. NF_NOERR) write (,) 'err in closing file',status
	623
	624
[524]	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