Context Navigation

readsulfate.F @ 1233

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

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"
	38	#include "temps.h"
[1146]	39	#include "clesphys.h"
	40	#include "iniprint.h"
[524]	41	c
	42	c Input:
	43	c ------
[782]	44	REAL r_day ! Day of integration
[524]	45	LOGICAL first ! First timestep
	46	! (and therefore initialization necessary)
	47	c
	48	c Output:
	49	c -------
[782]	50	REAL sulfate_p(klon_omp,klev)
	51	REAL sulfate (klon, klev) ! Mass of sulfate (monthly mean data,
[524]	52	! from file) [ug SO4/m3]
	53	c
	54	c Local Variables:
	55	c ----------------
	56	INTEGER i, ig, k, it
[640]	57	INTEGER j, iday, ny, iyr, iyr1, iyr2
[524]	58	parameter (ny=jjm+1)
	59
[640]	60	CJLD INTEGER idec1, idec2 ! The two decadal data read ini
[524]	61	CHARACTER*4 cyear
	62
	63	INTEGER im, day1, day2, im2
[782]	64	REAL so4_1(iim, jjm+1, klev, 12)
	65	REAL so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions
[524]	66
[782]	67	REAL, allocatable,save :: so4(:, :, :) ! SO4 in right dimension
	68	REAL, allocatable,save :: so4_out(:, :)
[766]	69	c$OMP THREADPRIVATE(so4,so4_out)
[524]	70
	71	LOGICAL lnewday
	72	LOGICAL lonlyone
	73	PARAMETER (lonlyone=.FALSE.)
[766]	74	logical,save :: first2=.true.
	75	c$OMP THREADPRIVATE(first2)
[524]	76
[766]	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
[776]	86	if (is_mpi_root) then
[1146]	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
[766]	96
[524]	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
[640]	129
[1146]	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)
[640]	137	ELSE
[1146]	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
[640]	147
	148	! Read in data:
[1146]	149	! a) from actual 10-yr-period
[524]	150
[1146]	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)
[524]	165
	166
	167	! If to read two decades:
[1146]	168	IF (.NOT.lonlyone) THEN
[524]	169
	170	! b) from the next following one
[1146]	171	WRITE(cyear,'(I4)') iyr2
	172	WRITE(,) 'getso4 iyr=', iyr,' ',cyear
	173	CALL getso4fromfile(cyear, so4_2)
[640]	174
	175
[524]	176	! Interpolate linarily to the actual year:
[1146]	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)
[640]	182	. - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
[524]	183	. * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
[1146]	184	ENDDO
	185	ENDDO
	186	ENDDO
	187	ENDDO
	188
	189
	190	ENDIF !lonlyone
	191	ENDIF !aer_type
[524]	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'
[804]	214	so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	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
[640]	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.
[524]	281	DO k=1,klev
	282	DO i=1,klon
[640]	283	CJLD sulfate(i,k) = sulfate(i,k)*masse_so4
	284	CJLD . /n_avogadro*1.e+12
[524]	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?)
[766]	301
	302	endif ! phy_rank==0
[524]	303
[766]	304	c$OMP END MASTER
[782]	305	call Scatter(sulfate,sulfate_p)
[766]	306
[524]	307	RETURN
	308	END
	309
	310
	311
	312
	313
	314	c-----------------------------------------------------------------------------
	315	c Read in /calculate pre-industrial values of sulfate
	316	c-----------------------------------------------------------------------------
	317
[766]	318	SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p)
[776]	319	USE dimphy, ONLY : klev
	320	USE mod_grid_phy_lmdz, klon=>klon_glo
	321	USE mod_phys_lmdz_para
[524]	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 ------
[782]	352	REAL r_day ! Day of integration
[524]	353	LOGICAL first ! First timestep
	354	! (and therefore initialization necessary)
	355	c
	356	c Output:
	357	c -------
[782]	358	REAL pi_sulfate_p (klon_omp, klev)
[766]	359
[782]	360	REAL pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data,
[766]	361	! from fil
[524]	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
[782]	369	INTEGER im, day1, day2, im2
	370	REAL pi_so4_1(iim, jjm+1, klev, 12)
[766]	371
[782]	372	REAL, allocatable,save :: pi_so4(:, :, :) ! SO4 in right dimension
	373	REAL, allocatable,save :: pi_so4_out(:, :)
[766]	374	c$OMP THREADPRIVATE(pi_so4,pi_so4_out)
[524]	375
	376	CHARACTER*4 cyear
	377	LOGICAL lnewday
[766]	378	logical,save :: first2=.true.
	379	c$OMP THREADPRIVATE(first2)
[524]	380
[766]	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
[776]	390	if (is_mpi_root) then
[766]	391
[524]	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'
[804]	452	pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it)
[524]	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
[640]	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]
[524]	518	DO k=1,klev
	519	DO i=1,klon
[640]	520	CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
	521	CJLD . /n_avogadro*1.e+12
[524]	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?
[766]	535
[776]	536	endif ! is_mpi_root==0
[766]	537
	538	c$OMP END MASTER
[782]	539	call Scatter(pi_sulfate,pi_sulfate_p)
[766]	540
[524]	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)
[1146]	559	USE dimphy
[524]	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
[782]	572	REAL so4mth(iim, ny, klev)
	573	REAL so4(iim, ny, klev, 12)
[524]	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
[782]	617
[524]	618	IF (STATUS .NE. NF_NOERR) write (,) 'err in read ',status
[782]	619
	620	#ifdef NC_DOUBLE
[801]	621	status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth)
[782]	622	#else
[801]	623	status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth)
[782]	624	#endif
[524]	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)
[782]	641	IF (STATUS .NE. NF_NOERR) write (,) 'err in closing file',status
	642
	643
[524]	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