Context Navigation

inifis.F @ 775

Last change on this file since 775 was 728, checked in by jleconte, 13 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

File size: 19.8 KB

Rev	Line
[135]	1	SUBROUTINE inifis(ngrid,nlayer,
	2	$ day_ini,pdaysec,ptimestep,
	3	$ plat,plon,parea,
	4	$ prad,pg,pr,pcpp)
[253]	5
	6	use radinc_h, only : naerkind
[374]	7	use datafile_mod, only: datadir
[253]	8
[135]	9	!=======================================================================
	10	!
	11	! purpose:
	12	! -------
	13	!
	14	! Initialisation for the physical parametrisations of the LMD
[305]	15	! Generic Model.
[135]	16	!
	17	! author: Frederic Hourdin 15 / 10 /93
	18	! -------
	19	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
	20	! Ehouarn Millour (oct. 2008) tracers are now identified
	21	! by their names and may not be contiguously
	22	! stored in the q(:,:,:,:) array
	23	! E.M. (june 2009) use getin routine to load parameters
	24	!
	25	!
	26	! arguments:
	27	! ----------
	28	!
	29	! input:
	30	! ------
	31	!
	32	! ngrid Size of the horizontal grid.
	33	! All internal loops are performed on that grid.
	34	! nlayer Number of vertical layers.
	35	! pdayref Day of reference for the simulation
	36	! pday Number of days counted from the North. Spring
	37	! equinoxe.
	38	!
	39	!=======================================================================
	40	!
	41	!-----------------------------------------------------------------------
	42	! declarations:
	43	! -------------
[374]	44	use datafile_mod, only: datadir
[135]	45	! to use 'getin'
	46	USE ioipsl_getincom
	47	IMPLICIT NONE
	48	#include "dimensions.h"
	49	#include "dimphys.h"
	50	#include "planete.h"
	51	#include "comcstfi.h"
	52	#include "comsaison.h"
	53	#include "comdiurn.h"
	54	#include "comgeomfi.h"
	55	#include "callkeys.h"
	56	#include "surfdat.h"
	57
	58
[305]	59
[135]	60	REAL prad,pg,pr,pcpp,pdaysec,ptimestep
	61
	62	INTEGER ngrid,nlayer
	63	REAL plat(ngrid),plon(ngrid),parea(ngridmx)
	64	integer day_ini
	65	INTEGER ig,ierr
	66
	67	EXTERNAL iniorbit,orbite
	68	EXTERNAL SSUM
	69	REAL SSUM
	70
	71	CHARACTER ch1*12
	72	CHARACTER ch80*80
	73
	74	logical chem, h2o
	75	logical :: parameter, doubleq=.false.
	76
[253]	77	real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
	78
[135]	79	rad=prad
	80	daysec=pdaysec
	81	dtphys=ptimestep
	82	cpp=pcpp
	83	g=pg
	84	r=pr
	85	rcp=r/cpp
	86
	87	avocado = 6.02214179e23 ! added by RW
	88
	89	! --------------------------------------------------------
	90	! The usual Tests
	91	! --------------
	92	IF (nlayer.NE.nlayermx) THEN
	93	PRINT*,'STOP in inifis'
	94	PRINT*,'Probleme de dimensions :'
	95	PRINT*,'nlayer = ',nlayer
	96	PRINT*,'nlayermx = ',nlayermx
	97	STOP
	98	ENDIF
	99
	100	IF (ngrid.NE.ngridmx) THEN
	101	PRINT*,'STOP in inifis'
	102	PRINT*,'Probleme de dimensions :'
	103	PRINT*,'ngrid = ',ngrid
	104	PRINT*,'ngridmx = ',ngridmx
	105	STOP
	106	ENDIF
	107
	108	! --------------------------------------------------------------
	109	! Reading the "callphys.def" file controlling some key options
	110	! --------------------------------------------------------------
	111
	112	! check that 'callphys.def' file is around
	113	OPEN(99,file='callphys.def',status='old',form='formatted'
	114	& ,iostat=ierr)
	115	CLOSE(99)
	116
	117	IF(ierr.EQ.0) THEN
	118	PRINT*
	119	PRINT*
	120	PRINT*,'--------------------------------------------'
	121	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
	122	PRINT*,'--------------------------------------------'
	123
[374]	124	write(,) "Directory where external input files are:"
	125	! default 'datadir' is set in "datadir_mod"
	126	call getin("datadir",datadir) ! default path
	127	write(,) " datadir = ",trim(datadir)
	128
[135]	129	write(,) "Run with or without tracer transport ?"
	130	tracer=.false. ! default value
	131	call getin("tracer",tracer)
	132	write(,) " tracer = ",tracer
	133
[728]	134	write(,) "Run with or without atm mass update ",
	135	& " due to tracer evaporation/condensation?"
	136	mass_redistrib=.false. ! default value
	137	call getin("mass_redistrib",mass_redistrib)
	138	write(,) " mass_redistrib = ",mass_redistrib
	139
[135]	140	write(,) "Diurnal cycle ?"
	141	write(,) "(if diurnal=false, diurnal averaged solar heating)"
	142	diurnal=.true. ! default value
	143	call getin("diurnal",diurnal)
	144	write(,) " diurnal = ",diurnal
	145
	146	write(,) "Seasonal cycle ?"
	147	write(,) "(if season=false, Ls stays constant, to value ",
	148	& "set in 'start'"
	149	season=.true. ! default value
	150	call getin("season",season)
	151	write(,) " season = ",season
	152
	153	write(,) "Tidally resonant rotation ?"
	154	tlocked=.false. ! default value
	155	call getin("tlocked",tlocked)
	156	write(,) "tlocked = ",tlocked
	157
	158	! Test of incompatibility:
	159	! if tlocked, then diurnal should be false
	160	if (tlocked.and.diurnal) then
	161	print*,'If diurnal=true, we should turn off tlocked.'
	162	stop
	163	endif
	164
	165	write(,) "Tidal resonance ratio ?"
	166	nres=0 ! default value
	167	call getin("nres",nres)
	168	write(,) "nres = ",nres
	169
	170	write(,) "Write some extra output to the screen ?"
	171	lwrite=.false. ! default value
	172	call getin("lwrite",lwrite)
	173	write(,) " lwrite = ",lwrite
	174
	175	write(,) "Save statistics in file stats.nc ?"
	176	callstats=.true. ! default value
	177	call getin("callstats",callstats)
	178	write(,) " callstats = ",callstats
	179
	180	write(,) "Test energy conservation of model physics ?"
	181	enertest=.false. ! default value
	182	call getin("enertest",enertest)
	183	write(,) " enertest = ",enertest
	184
[538]	185	write(,) "Check to see if cpp values used match gases.def ?"
	186	check_cpp_match=.true. ! default value
	187	call getin("check_cpp_match",check_cpp_match)
	188	write(,) " check_cpp_match = ",check_cpp_match
	189
[135]	190	write(,) "call radiative transfer ?"
	191	callrad=.true. ! default value
	192	call getin("callrad",callrad)
	193	write(,) " callrad = ",callrad
	194
	195	write(,) "call correlated-k radiative transfer ?"
	196	corrk=.true. ! default value
	197	call getin("corrk",corrk)
	198	write(,) " corrk = ",corrk
	199
	200	write(,) "call gaseous absorption in the visible bands?",
	201	& "(matters only if callrad=T)"
	202	callgasvis=.false. ! default value
	203	call getin("callgasvis",callgasvis)
	204	write(,) " callgasvis = ",callgasvis
[538]	205
[526]	206	write(,) "call continuum opacities in radiative transfer ?",
	207	& "(matters only if callrad=T)"
[538]	208	Continuum=.true. ! default value
[526]	209	call getin("Continuum",Continuum)
	210	write(,) " Continuum = ",Continuum
[538]	211
[716]	212	write(,) "use analytic function for H2O continuum ?"
	213	H2Ocont_simple=.false. ! default value
	214	call getin("H2Ocont_simple",H2Ocont_simple)
	215	write(,) " H2Ocont_simple = ",H2Ocont_simple
[538]	216
[135]	217	write(,) "call turbulent vertical diffusion ?"
	218	calldifv=.true. ! default value
	219	call getin("calldifv",calldifv)
	220	write(,) " calldifv = ",calldifv
	221
[596]	222	write(,) "use turbdiff instead of vdifc ?"
	223	UseTurbDiff=.true. ! default value
	224	call getin("UseTurbDiff",UseTurbDiff)
	225	write(,) " UseTurbDiff = ",UseTurbDiff
	226
[135]	227	write(,) "call convective adjustment ?"
	228	calladj=.true. ! default value
	229	call getin("calladj",calladj)
	230	write(,) " calladj = ",calladj
	231
	232	write(,) "call CO2 condensation ?"
[716]	233	co2cond=.false. ! default value
[135]	234	call getin("co2cond",co2cond)
	235	write(,) " co2cond = ",co2cond
[650]	236	! Test of incompatibility
	237	if (co2cond.and.(.not.tracer)) then
	238	print*,'We need a CO2 ice tracer to condense CO2'
	239	call abort
[716]	240	endif
	241
[253]	242	write(,) "CO2 supersaturation level ?"
	243	co2supsat=1.0 ! default value
	244	call getin("co2supsat",co2supsat)
	245	write(,) " co2supsat = ",co2supsat
	246
	247	write(,) "Radiative timescale for Newtonian cooling ?"
	248	tau_relax=30. ! default value
	249	call getin("tau_relax",tau_relax)
	250	write(,) " tau_relax = ",tau_relax
	251	tau_relax=tau_relax243600 ! convert Earth days --> seconds
	252
[135]	253	write(,)"call thermal conduction in the soil ?"
	254	callsoil=.true. ! default value
	255	call getin("callsoil",callsoil)
	256	write(,) " callsoil = ",callsoil
	257
[253]	258	write(,)"Rad transfer is computed every iradia",
[135]	259	& " physical timestep"
	260	iradia=1 ! default value
	261	call getin("iradia",iradia)
	262	write(,)" iradia = ",iradia
	263
	264	write(,)"Rayleigh scattering ?"
	265	rayleigh=.false.
	266	call getin("rayleigh",rayleigh)
	267	write(,)" rayleigh = ",rayleigh
	268
[253]	269	write(,) "Use blackbody for stellar spectrum ?"
	270	stelbbody=.false. ! default value
	271	call getin("stelbbody",stelbbody)
	272	write(,) " stelbbody = ",stelbbody
[135]	273
[253]	274	write(,) "Stellar blackbody temperature ?"
	275	stelTbb=5800.0 ! default value
	276	call getin("stelTbb",stelTbb)
	277	write(,) " stelTbb = ",stelTbb
	278
[135]	279	write(,)"Output mean OLR in 1D?"
	280	meanOLR=.false.
	281	call getin("meanOLR",meanOLR)
	282	write(,)" meanOLR = ",meanOLR
	283
	284	write(,)"Output spectral OLR in 3D?"
	285	specOLR=.false.
	286	call getin("specOLR",specOLR)
	287	write(,)" specOLR = ",specOLR
	288
[253]	289	write(,)"Operate in kastprof mode?"
	290	kastprof=.false.
	291	call getin("kastprof",kastprof)
	292	write(,)" kastprof = ",kastprof
	293
[596]	294	write(,)"Uniform absorption in radiative transfer?"
[305]	295	graybody=.false.
	296	call getin("graybody",graybody)
	297	write(,)" graybody = ",graybody
[253]	298
	299	! Test of incompatibility:
	300	! if kastprof used, we must be in 1D
	301	if (kastprof.and.(ngridmx.gt.1)) then
	302	print*,'kastprof can only be used in 1D!'
	303	call abort
	304	endif
	305
	306	write(,)"Stratospheric temperature for kastprof mode?"
	307	Tstrat=167.0
	308	call getin("Tstrat",Tstrat)
	309	write(,)" Tstrat = ",Tstrat
	310
	311	write(,)"Remove lower boundary?"
	312	noradsurf=.false.
	313	call getin("noradsurf",noradsurf)
	314	write(,)" noradsurf = ",noradsurf
	315
	316	! Tests of incompatibility:
	317	if (noradsurf.and.callsoil) then
	318	print*,'noradsurf not compatible with soil scheme!'
	319	call abort
	320	endif
[589]	321	!if (noradsurf.and.calldifv) then
	322	! print*,'noradsurf not compatible with a boundary layer!'
	323	! call abort
	324	!endif
[253]	325
	326	write(,)"Use Newtonian cooling for radiative transfer?"
	327	newtonian=.false.
	328	call getin("newtonian",newtonian)
	329	write(,)" newtonian = ",newtonian
	330
	331	! Tests of incompatibility:
	332	if (newtonian.and.corrk) then
	333	print*,'newtonian not compatible with correlated-k!'
	334	call abort
	335	endif
	336	if (newtonian.and.calladj) then
	337	print*,'newtonian not compatible with adjustment!'
	338	call abort
	339	endif
	340	if (newtonian.and.calldifv) then
	341	print*,'newtonian not compatible with a boundary layer!'
	342	call abort
	343	endif
	344
	345	write(,)"Test physics timescale in 1D?"
	346	testradtimes=.false.
	347	call getin("testradtimes",testradtimes)
	348	write(,)" testradtimes = ",testradtimes
	349
	350	! Test of incompatibility:
	351	! if testradtimes used, we must be in 1D
	352	if (testradtimes.and.(ngridmx.gt.1)) then
	353	print*,'testradtimes can only be used in 1D!'
	354	call abort
	355	endif
	356
[135]	357	write(,)"Default planetary temperature?"
	358	tplanet=215.0
	359	call getin("tplanet",tplanet)
	360	write(,)" tplanet = ",tplanet
	361
	362	write(,)"Which star?"
	363	startype=1 ! default value = Sol
	364	call getin("startype",startype)
	365	write(,)" startype = ",startype
	366
	367	write(,)"Value of stellar flux at 1 AU?"
	368	Fat1AU=1356.0 ! default value = Sol today
	369	call getin("Fat1AU",Fat1AU)
	370	write(,)" Fat1AU = ",Fat1AU
	371
	372
	373	! TRACERS:
	374
[253]	375	write(,)"Varying H2O cloud fraction?"
	376	CLFvarying=.false. ! default value
	377	call getin("CLFvarying",CLFvarying)
	378	write(,)" CLFvarying = ",CLFvarying
	379
	380	write(,)"Value of fixed H2O cloud fraction?"
	381	CLFfixval=1.0 ! default value
	382	call getin("CLFfixval",CLFfixval)
	383	write(,)" CLFfixval = ",CLFfixval
	384
[728]	385	write(,)"fixed radii for Cloud particles?"
	386	radfixed=.false. ! default value
	387	call getin("radfixed",radfixed)
	388	write(,)" radfixed = ",radfixed
	389
	390	if(kastprof)then
	391	radfixed=.true.
	392	endif
	393
	394	write(,)"Number mixing ratio of CO2 ice particles:"
	395	Nmix_co2=1.e6 ! default value
[135]	396	call getin("Nmix_co2",Nmix_co2)
	397	write(,)" Nmix_co2 = ",Nmix_co2
	398
[726]	399	! write(,)"Number of radiatively active aerosols:"
	400	! naerkind=0. ! default value
	401	! call getin("naerkind",naerkind)
	402	! write(,)" naerkind = ",naerkind
	403
[253]	404	write(,)"Opacity of dust (if used):"
	405	dusttau=0. ! default value
	406	call getin("dusttau",dusttau)
	407	write(,)" dusttau = ",dusttau
	408
[726]	409	write(,)"Radiatively active CO2 aerosols?"
	410	aeroco2=.false. ! default value
	411	call getin("aeroco2",aeroco2)
	412	write(,)" aeroco2 = ",aeroco2
[253]	413
[726]	414	write(,)"Fixed CO2 aerosol distribution?"
	415	aerofixco2=.false. ! default value
	416	call getin("aerofixco2",aerofixco2)
	417	write(,)" aerofixco2 = ",aerofixco2
	418
	419	write(,)"Radiatively active water ice?"
	420	aeroh2o=.false. ! default value
	421	call getin("aeroh2o",aeroh2o)
	422	write(,)" aeroh2o = ",aeroh2o
	423
	424	write(,)"Fixed H2O aerosol distribution?"
	425	aerofixh2o=.false. ! default value
	426	call getin("aerofixh2o",aerofixh2o)
	427	write(,)" aerofixh2o = ",aerofixh2o
	428
	429	write(,)"Radiatively active sulfuric acid aersols?"
	430	aeroh2so4=.false. ! default value
	431	call getin("aeroh2so4",aeroh2so4)
	432	write(,)" aeroh2so4 = ",aeroh2so4
	433
[305]	434	write(,)"Cloud pressure level (with kastprof only):"
	435	cloudlvl=0. ! default value
	436	call getin("cloudlvl",cloudlvl)
	437	write(,)" cloudlvl = ",cloudlvl
	438
[135]	439	write(,)"Is the variable gas species radiatively active?"
[305]	440	Tstrat=167.0
[135]	441	varactive=.false.
	442	call getin("varactive",varactive)
	443	write(,)" varactive = ",varactive
	444
	445	write(,)"Is the variable gas species distribution set?"
	446	varfixed=.false.
	447	call getin("varfixed",varfixed)
	448	write(,)" varfixed = ",varfixed
	449
	450	write(,)"What is the saturation % of the variable species?"
	451	satval=0.8
	452	call getin("satval",satval)
	453	write(,)" satval = ",satval
	454
	455	! Test of incompatibility:
[726]	456	! if no tracers, then the variable aerosols should be turned off
	457	if ((.not.tracer).and.(.not.aerofixco2)) then
	458	print*,'if tracers are off, you must turn aerofixco2 on.'
	459	print*, '(i.e., you must fix the aerosol distribution)'
	460	print*, 'This option was formerly known as aerofixed'
[135]	461	stop
	462	endif
[726]	463	if ((.not.tracer).and.(.not.aerofixh2o)) then
	464	print*,'if tracers are off, you must turn aerofixh2o on.'
	465	print*, '(i.e., you must fix the aerosol distribution)'
	466	print*, 'This option was formerly known as aerofixed'
	467	stop
	468	endif
[135]	469
[726]	470
[135]	471	! Test of incompatibility:
	472	! if varactive, then varfixed should be false
	473	if (varactive.and.varfixed) then
	474	print*,'if varactive, varfixed must be OFF!'
	475	stop
	476	endif
	477
	478	write(,) "Gravitationnal sedimentation ?"
[716]	479	sedimentation=.false. ! default value
[135]	480	call getin("sedimentation",sedimentation)
	481	write(,) " sedimentation = ",sedimentation
	482
	483	write(,) "Compute water cycle ?"
	484	water=.false. ! default value
	485	call getin("water",water)
	486	write(,) " water = ",water
	487
[622]	488	! Test of incompatibility:
	489	! if water is true, there should be at least a tracer
	490	if (water.and.(.not.tracer)) then
	491	print*,'if water is ON, tracer must be ON too!'
	492	stop
	493	endif
	494
[135]	495	write(,) "Include water condensation ?"
	496	watercond=.false. ! default value
	497	call getin("watercond",watercond)
	498	write(,) " watercond = ",watercond
	499
[622]	500	! Test of incompatibility:
	501	! if watercond is used, then water should be used too
	502	if (watercond.and.(.not.water)) then
	503	print*,'if watercond is used, water should be used too'
	504	stop
	505	endif
	506
[135]	507	write(,) "Include water precipitation ?"
	508	waterrain=.false. ! default value
	509	call getin("waterrain",waterrain)
	510	write(,) " waterrain = ",waterrain
	511
[253]	512	write(,) "Include surface hydrology ?"
	513	hydrology=.false. ! default value
	514	call getin("hydrology",hydrology)
	515	write(,) " hydrology = ",hydrology
[135]	516
[253]	517	write(,) "Evolve surface water sources ?"
	518	sourceevol=.false. ! default value
	519	call getin("sourceevol",sourceevol)
	520	write(,) " sourceevol = ",sourceevol
[135]	521
[486]	522	write(,) "Ice evolution timestep ?"
	523	icetstep=100.0 ! default value
	524	call getin("icetstep",icetstep)
	525	write(,) " icetstep = ",icetstep
	526
[253]	527	write(,) "Snow albedo ?"
	528	albedosnow=0.5 ! default value
	529	call getin("albedosnow",albedosnow)
	530	write(,) " albedosnow = ",albedosnow
[135]	531
[253]	532	write(,) "Maximum ice thickness ?"
	533	maxicethick=2.0 ! default value
	534	call getin("maxicethick",maxicethick)
	535	write(,) " maxicethick = ",maxicethick
[135]	536
[253]	537	write(,) "Freezing point of seawater ?"
	538	Tsaldiff=-1.8 ! default value
	539	call getin("Tsaldiff",Tsaldiff)
	540	write(,) " Tsaldiff = ",Tsaldiff
[135]	541
[589]	542	write(,) "Does user want to force cpp and mugaz?"
	543	force_cpp=.false. ! default value
	544	call getin("force_cpp",force_cpp)
	545	write(,) " force_cpp = ",force_cpp
	546
	547	if (force_cpp) then
	548	mugaz = -99999.
	549	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
	550	call getin("mugaz",mugaz)
	551	IF (mugaz.eq.-99999.) THEN
	552	PRINT *, "mugaz must be set if force_cpp = T"
	553	STOP
	554	ELSE
	555	write(,) "mugaz=",mugaz
	556	ENDIF
	557	cpp = -99999.
	558	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
	559	call getin("cpp",cpp)
	560	IF (cpp.eq.-99999.) THEN
	561	PRINT *, "cpp must be set if force_cpp = T"
	562	STOP
	563	ELSE
	564	write(,) "cpp=",cpp
	565	ENDIF
	566	else
	567	mugaz=8.314*1000./pr
	568	endif
[253]	569	call su_gases
	570	call calc_cpp_mugaz
[135]	571
	572	PRINT*,'--------------------------------------------'
	573	PRINT*
	574	PRINT*
	575	ELSE
	576	write(,)
	577	write(,) 'Cannot read file callphys.def. Is it here ?'
	578	stop
	579	ENDIF
	580
	581	8000 FORMAT(t5,a12,l8)
	582	8001 FORMAT(t5,a12,i8)
	583
	584	PRINT*
	585	PRINT*,'inifis: daysec',daysec
	586	PRINT*
	587	PRINT*,'inifis: The radiative transfer is computed:'
	588	PRINT*,' each ',iradia,' physical time-step'
	589	PRINT,' or each ',iradiadtphys,' seconds'
	590	PRINT*
	591
	592
	593	!-----------------------------------------------------------------------
	594	! Some more initialization:
	595	! ------------------------
	596
	597	CALL SCOPY(ngrid,plon,1,long,1)
	598	CALL SCOPY(ngrid,plat,1,lati,1)
	599	CALL SCOPY(ngrid,parea,1,area,1)
	600	totarea=SSUM(ngridmx,area,1)
	601
	602	DO ig=1,ngrid
	603	sinlat(ig)=sin(plat(ig))
	604	coslat(ig)=cos(plat(ig))
	605	sinlon(ig)=sin(plon(ig))
	606	coslon(ig)=cos(plon(ig))
	607	ENDDO
	608
	609	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
	610
	611	RETURN
	612	END

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

Download in other formats:

Original Format