Context Navigation

inifis.F @ 1088

Last change on this file since 1088 was 1088, checked in by tnavarro, 11 years ago
Added a freedust mode, to be used with doubleq (for data assimilation). Dust is not lifted, not rescaled, dust opacity is predicted instead of being forced.
File size: 27.9 KB

Line
1	SUBROUTINE inifis(
2	$ ngrid,nlayer,nq
3	$ ,day_ini,pdaysec,ptimestep
4	$ ,plat,plon,parea
5	$ ,prad,pg,pr,pcpp
6	#ifdef MESOSCALE
7	#include "meso_inc/meso_inc_inifisinvar.F"
8	#endif
9	$ )
10	!
11	!=======================================================================
12	!
13	! purpose:
14	! -------
15	!
16	! Initialisation for the physical parametrisations of the LMD
17	! martian atmospheric general circulation modele.
18	!
19	! author: Frederic Hourdin 15 / 10 /93
20	! -------
21	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
22	! Ehouarn Millour (oct. 2008) tracers are now identified
23	! by their names and may not be contiguously
24	! stored in the q(:,:,:,:) array
25	! E.M. (june 2009) use getin routine to load parameters
26	! adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
27	!
28	!
29	! arguments:
30	! ----------
31	!
32	! input:
33	! ------
34	!
35	! ngrid Size of the horizontal grid.
36	! All internal loops are performed on that grid.
37	! nlayer Number of vertical layers.
38	! pdayref Day of reference for the simulation
39	! pday Number of days counted from the North. Spring
40	! equinoxe.
41	!
42	!=======================================================================
43	!
44	!-----------------------------------------------------------------------
45	! declarations:
46	! -------------
47	! to use 'getin'
48	USE ioipsl_getincom, only : getin
49	use tracer_mod, only : nqmx, nuice_sed, ccn_factor
50	use comsoil_h, only: ini_comsoil_h
51	#ifdef MESOSCALE
52	use comsoil_h, only: volcapa !!MESOSCALE -- needed to fill volcapa
53	#endif
54	use comgeomfi_h, only: long, lati, area, totarea
55	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
56	use surfdat_h, only: ini_surfdat_h, albedo_h2o_ice, inert_h2o_ice,
57	& frost_albedo_threshold
58	use comsaison_h, only: ini_comsaison_h
59	use slope_mod, only: ini_slope_mod
60	use dimradmars_mod, only: ini_dimradmars_mod
61	use yomaer_h,only: ini_yomaer_h, tauvis
62	use yomlw_h, only: ini_yomlw_h
63	use conc_mod, only: ini_conc_mod
64	IMPLICIT NONE
65	#include "dimensions.h"
66	#include "dimphys.h"
67	#include "planete.h"
68	#include "comcstfi.h"
69	!#include "comsaison.h"
70	!#include "comdiurn.h"
71	!#include "comgeomfi.h"
72	#include "callkeys.h"
73	!#include "surfdat.h"
74	!#include "dimradmars.h"
75	!#include "yomaer.h"
76	#include "datafile.h"
77	!#include "slope.h"
78	#include "microphys.h"
79	!#include "tracer.h"
80	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
81	#include"scatterers.h"
82	#ifdef MESOSCALE
83	!#include "comsoil.h" !!MESOSCALE -- needed to fill volcapa
84	#include "meso_inc/meso_inc_inifisvar.F"
85	#endif
86	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec
87
88	REAL,INTENT(IN) :: ptimestep
89	INTEGER,INTENT(IN) :: day_ini
90
91	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
92	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
93	INTEGER ig,ierr
94
95	! EXTERNAL iniorbit,orbite
96	EXTERNAL SSUM
97	REAL SSUM
98
99	CHARACTER ch1*12
100	CHARACTER ch80*80
101
102	! logical chem, h2o
103
104	! chem = .false.
105	! h2o = .false.
106
107	rad=prad
108	cpp=pcpp
109	g=pg
110	r=pr
111	rcp=r/cpp
112	daysec=pdaysec
113	dtphys=ptimestep
114
115	nqmx=nq
116
117	#ifdef MESOSCALE
118	#include "meso_inc/meso_inc_inifisini.F"
119	#endif
120
121	! --------------------------------------------------------
122	! The usual Tests
123	! --------------
124	! IF (nlayer.NE.nlayermx) THEN
125	! PRINT*,'STOP in inifis'
126	! PRINT*,'Probleme de dimensions :'
127	! PRINT*,'nlayer = ',nlayer
128	! PRINT*,'nlayermx = ',nlayermx
129	! STOP
130	! ENDIF
131
132	! IF (ngrid.NE.ngridmx) THEN
133	! PRINT*,'STOP in inifis'
134	! PRINT*,'Probleme de dimensions :'
135	! PRINT*,'ngrid = ',ngrid
136	! PRINT*,'ngridmx = ',ngridmx
137	! STOP
138	! ENDIF
139
140	! --------------------------------------------------------------
141	! Reading the "callphys.def" file controlling some key options
142	! --------------------------------------------------------------
143
144	! check that 'callphys.def' file is around
145	OPEN(99,file='callphys.def',status='old',form='formatted'
146	& ,iostat=ierr)
147	CLOSE(99)
148
149	IF(ierr.EQ.0) THEN
150	PRINT*
151	PRINT*
152	PRINT*,'--------------------------------------------'
153	PRINT*,' inifis: Parameters for the physics (callphys.def)'
154	PRINT*,'--------------------------------------------'
155
156	write(,) "Directory where external input files are:"
157	datafile="/u/forget/WWW/datagcm/datafile"
158	call getin("datadir",datafile) ! default path
159	write(,) " datafile = ",trim(datafile)
160
161	write(,) "Run with or without tracer transport ?"
162	tracer=.false. ! default value
163	call getin("tracer",tracer)
164	write(,) " tracer = ",tracer
165
166	write(,) "Diurnal cycle ?"
167	write(,) "(if diurnal=False, diurnal averaged solar heating)"
168	diurnal=.true. ! default value
169	call getin("diurnal",diurnal)
170	write(,) " diurnal = ",diurnal
171
172	write(,) "Seasonal cycle ?"
173	write(,) "(if season=False, Ls stays constant, to value ",
174	& "set in 'start'"
175	season=.true. ! default value
176	call getin("season",season)
177	write(,) " season = ",season
178
179	write(,) "Write some extra output to the screen ?"
180	lwrite=.false. ! default value
181	call getin("lwrite",lwrite)
182	write(,) " lwrite = ",lwrite
183
184	write(,) "Save statistics in file stats.nc ?"
185	#ifdef MESOSCALE
186	callstats=.false. ! default value
187	#else
188	callstats=.true. ! default value
189	#endif
190	call getin("callstats",callstats)
191	write(,) " callstats = ",callstats
192
193	write(,) "Save EOF profiles in file 'profiles' for ",
194	& "Climate Database?"
195	calleofdump=.false. ! default value
196	call getin("calleofdump",calleofdump)
197	write(,) " calleofdump = ",calleofdump
198
199	write(,) "Dust scenario: 1=constant dust (read from startfi",
200	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
201	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
202	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
203	iaervar=3 ! default value
204	call getin("iaervar",iaervar)
205	write(,) " iaervar = ",iaervar
206
207	write(,) "Reference (visible) dust opacity at 610 Pa ",
208	& "(matters only if iaervar=1)"
209	! NB: default value of tauvis is set/read in startfi.nc file
210	call getin("tauvis",tauvis)
211	write(,) " tauvis = ",tauvis
212
213	write(,) "Dust vertical distribution:"
214	write(,) "(=1 top set by topdustref parameter;",
215	& " =2 Viking scenario; =3 MGS scenario)"
216	iddist=3 ! default value
217	call getin("iddist",iddist)
218	write(,) " iddist = ",iddist
219
220	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
221	topdustref= 90.0 ! default value
222	call getin("topdustref",topdustref)
223	write(,) " topdustref = ",topdustref
224
225	write(,) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
226	tke_heat_flux=0. ! default value
227	call getin("tke_heat_flux",tke_heat_flux)
228	write(,) " tke_heat_flux = ",tke_heat_flux
229	write(,) " 0 means the usual schemes are computing"
230
231	write(,) "call radiative transfer ?"
232	callrad=.true. ! default value
233	call getin("callrad",callrad)
234	write(,) " callrad = ",callrad
235
236	write(,) "call slope insolation scheme ?",
237	& "(matters only if callrad=T)"
238	#ifdef MESOSCALE
239	callslope=.true. ! default value
240	#else
241	callslope=.false. ! default value (not supported yet)
242	#endif
243	call getin("callslope",callslope)
244	write(,) " callslope = ",callslope
245
246	write(,) "call NLTE radiative schemes ?",
247	& "(matters only if callrad=T)"
248	callnlte=.false. ! default value
249	call getin("callnlte",callnlte)
250	write(,) " callnlte = ",callnlte
251
252	nltemodel=0 !default value
253	write(,) "NLTE model?"
254	write(,) "0 -> old model, static O"
255	write(,) "1 -> old model, dynamic O"
256	write(,) "2 -> new model"
257	write(,) "(matters only if callnlte=T)"
258	call getin("nltemodel",nltemodel)
259	write(,) " nltemodel = ",nltemodel
260
261	write(,) "call CO2 NIR absorption ?",
262	& "(matters only if callrad=T)"
263	callnirco2=.false. ! default value
264	call getin("callnirco2",callnirco2)
265	write(,) " callnirco2 = ",callnirco2
266
267	write(,) "New NIR NLTE correction ?",
268	$ "0-> old model (no correction)",
269	$ "1-> new correction",
270	$ "(matters only if callnirco2=T)"
271	#ifdef MESOSCALE
272	nircorr=0 !default value. this is OK below 60 km.
273	#else
274	nircorr=0 !default value
275	#endif
276	call getin("nircorr",nircorr)
277	write(,) " nircorr = ",nircorr
278
279	write(,) "call turbulent vertical diffusion ?"
280	calldifv=.true. ! default value
281	call getin("calldifv",calldifv)
282	write(,) " calldifv = ",calldifv
283
284	write(,) "call thermals ?"
285	calltherm=.false. ! default value
286	call getin("calltherm",calltherm)
287	write(,) " calltherm = ",calltherm
288
289	write(,) "call convective adjustment ?"
290	calladj=.true. ! default value
291	call getin("calladj",calladj)
292	write(,) " calladj = ",calladj
293
294	if (calltherm .and. (.not. calladj)) then
295	print*,'Convadj has to be activated when using thermals'
296	stop
297	endif
298
299	write(,) "call Richardson-based surface layer ?"
300	callrichsl=.false. ! default value
301	call getin("callrichsl",callrichsl)
302	write(,) " callrichsl = ",callrichsl
303
304	if (calltherm .and. .not.callrichsl) then
305	print*,'WARNING WARNING WARNING'
306	print*,'if calltherm=T we strongly advise that '
307	print*,'you use the new surface layer scheme '
308	print*,'by setting callrichsl=T '
309	endif
310
311	if (calladj .and. callrichsl .and. (.not. calltherm)) then
312	print*,'You should not be calling the convective adjustment
313	& scheme with the Richardson surface-layer and without the thermals
314	&. This approach is not
315	& physically consistent and can lead to unrealistic friction
316	& values.'
317	print*,'If you want to use the Ri. surface-layer, either
318	& activate thermals OR de-activate the convective adjustment.'
319	stop
320	endif
321
322	write(,) "call CO2 condensation ?"
323	callcond=.true. ! default value
324	call getin("callcond",callcond)
325	write(,) " callcond = ",callcond
326
327	write(,)"call thermal conduction in the soil ?"
328	callsoil=.true. ! default value
329	call getin("callsoil",callsoil)
330	write(,) " callsoil = ",callsoil
331
332
333	write(,)"call Lott's gravity wave/subgrid topography ",
334	& "scheme ?"
335	calllott=.true. ! default value
336	call getin("calllott",calllott)
337	write(,)" calllott = ",calllott
338
339
340	write(,)"rad.transfer is computed every iradia",
341	& " physical timestep"
342	iradia=1 ! default value
343	call getin("iradia",iradia)
344	write(,)" iradia = ",iradia
345
346
347	write(,)"Output of the exchange coefficient mattrix ?",
348	& "(for diagnostics only)"
349	callg2d=.false. ! default value
350	call getin("callg2d",callg2d)
351	write(,)" callg2d = ",callg2d
352
353	write(,)"Rayleigh scattering : (should be .false. for now)"
354	rayleigh=.false.
355	call getin("rayleigh",rayleigh)
356	write(,)" rayleigh = ",rayleigh
357
358
359	! TRACERS:
360
361	! dustbin
362	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
363	dustbin=0 ! default value
364	call getin("dustbin",dustbin)
365	write(,)" dustbin = ",dustbin
366	! active
367	write(,)"Radiatively active dust ? (matters if dustbin>0)"
368	active=.false. ! default value
369	call getin("active",active)
370	write(,)" active = ",active
371
372	! Test of incompatibility:
373	! if active is used, then dustbin should be > 0
374
375	if (active.and.(dustbin.lt.1)) then
376	print*,'if active is used, then dustbin should > 0'
377	stop
378	endif
379	! doubleq
380	write(,)"use mass and number mixing ratios to predict",
381	& " dust size ?"
382	doubleq=.false. ! default value
383	call getin("doubleq",doubleq)
384	write(,)" doubleq = ",doubleq
385	! submicron
386	submicron=.false. ! default value
387	call getin("submicron",submicron)
388	write(,)" submicron = ",submicron
389
390	! Test of incompatibility:
391	! if doubleq is used, then dustbin should be 2
392
393	if (doubleq.and.(dustbin.ne.2)) then
394	print*,'if doubleq is used, then dustbin should be 2'
395	stop
396	endif
397	if (doubleq.and.submicron.and.(nq.LT.3)) then
398	print*,'If doubleq is used with a submicron tracer,'
399	print*,' then the number of tracers has to be'
400	print*,' larger than 3.'
401	stop
402	endif
403
404	! lifting
405	write(,)"dust lifted by GCM surface winds ?"
406	lifting=.false. ! default value
407	call getin("lifting",lifting)
408	write(,)" lifting = ",lifting
409
410	! Test of incompatibility:
411	! if lifting is used, then dustbin should be > 0
412
413	if (lifting.and.(dustbin.lt.1)) then
414	print*,'if lifting is used, then dustbin should > 0'
415	stop
416	endif
417
418	! free evolving dust
419	! freedust=true just says that there is no lifting and no dust opacity scaling.
420	write(,)"dust lifted by GCM surface winds ?"
421	freedust=.false. ! default value
422	call getin("freedust",freedust)
423	write(,)" freedust = ",freedust
424	if (freedust.and..not.doubleq) then
425	print*,'freedust should be used with doubleq !'
426	stop
427	endif
428	if (freedust.and.lifting) then
429	print*,'if freedust is used, then lifting should not be used'
430	print*,'lifting forced to false !!'
431	lifting=.false.
432	endif
433
434	! callddevil
435	write(,)" dust lifted by dust devils ?"
436	callddevil=.false. !default value
437	call getin("callddevil",callddevil)
438	write(,)" callddevil = ",callddevil
439
440	! Test of incompatibility:
441	! if dustdevil is used, then dustbin should be > 0
442
443	if (callddevil.and.(dustbin.lt.1)) then
444	print*,'if dustdevil is used, then dustbin should > 0'
445	stop
446	endif
447	! sedimentation
448	write(,) "Gravitationnal sedimentation ?"
449	sedimentation=.true. ! default value
450	call getin("sedimentation",sedimentation)
451	write(,) " sedimentation = ",sedimentation
452	! activice
453	write(,) "Radiatively active transported atmospheric ",
454	& "water ice ?"
455	activice=.false. ! default value
456	call getin("activice",activice)
457	write(,) " activice = ",activice
458	! water
459	write(,) "Compute water cycle ?"
460	water=.false. ! default value
461	call getin("water",water)
462	write(,) " water = ",water
463
464	! thermal inertia feedback
465	write(,) "Activate the thermal inertia feedback ?"
466	tifeedback=.false. ! default value
467	call getin("tifeedback",tifeedback)
468	write(,) " tifeedback = ",tifeedback
469
470	! Test of incompatibility:
471
472	if (tifeedback.and..not.water) then
473	print*,'if tifeedback is used,'
474	print*,'water should be used too'
475	stop
476	endif
477
478	if (tifeedback.and..not.callsoil) then
479	print*,'if tifeedback is used,'
480	print*,'callsoil should be used too'
481	stop
482	endif
483
484	if (activice.and..not.water) then
485	print*,'if activice is used, water should be used too'
486	stop
487	endif
488
489	if (water.and..not.tracer) then
490	print*,'if water is used, tracer should be used too'
491	stop
492	endif
493
494	! water ice clouds effective variance distribution for sedimentaion
495	write(,) "effective variance for water ice clouds ?"
496	nuice_sed=0.45
497	call getin("nuice_sed",nuice_sed)
498	write(,) "water_param nueff Sedimentation:", nuice_sed
499
500	! ccn factor if no scavenging
501	write(,) "water param CCN reduc. factor ?", ccn_factor
502	ccn_factor = 4.5
503	call getin("ccn_factor",ccn_factor)
504	write(,)" ccn_factor = ",ccn_factor
505	write(,)"Careful: only used when microphys=F, otherwise"
506	write(,)"the contact parameter is used instead;"
507
508	! microphys
509	write(,)"Microphysical scheme for water-ice clouds?"
510	microphys=.false. ! default value
511	call getin("microphys",microphys)
512	write(,)" microphys = ",microphys
513	if (freedust.and.microphys) then
514	print*,'freedust can not be used with microphys !'
515	print*,'(although that could be improved...)'
516	stop
517	endif
518
519	! microphysical parameter contact
520	write(,) "water contact parameter ?"
521	mteta = 0.95
522	call getin("mteta",mteta)
523	write(,) "mteta = ", mteta
524
525	! scavenging
526	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
527	scavenging=.false. ! default value
528	call getin("scavenging",scavenging)
529	write(,)" scavenging = ",scavenging
530
531
532	! Test of incompatibility:
533	! if scavenging is used, then dustbin should be > 0
534
535	if ((microphys.and..not.doubleq).or.
536	& (microphys.and..not.water)) then
537	print*,'if microphys is used, then doubleq,'
538	print*,'and water must be used!'
539	stop
540	endif
541	if (microphys.and..not.scavenging) then
542	print*,''
543	print*,'----------------WARNING-----------------'
544	print*,'microphys is used without scavenging !!!'
545	print*,'----------------WARNING-----------------'
546	print*,''
547	endif
548
549	if ((scavenging.and..not.microphys).or.
550	& (scavenging.and.(dustbin.lt.1))) then
551	print*,'if scavenging is used, then microphys'
552	print*,'must be used!'
553	stop
554	endif
555
556	! Test of incompatibility:
557
558	write(,) "Permanent water caps at poles ?",
559	& " .true. is RECOMMENDED"
560	write(,) "(with .true., North cap is a source of water ",
561	& "and South pole is a cold trap)"
562	caps=.true. ! default value
563	call getin("caps",caps)
564	write(,) " caps = ",caps
565
566	! albedo_h2o_ice
567	write(,) "water ice albedo ?"
568	albedo_h2o_ice=0.45
569	call getin("albedo_h2o_ice",albedo_h2o_ice)
570	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
571	! inert_h2o_ice
572	write(,) "water ice thermal inertia ?"
573	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
574	call getin("inert_h2o_ice",inert_h2o_ice)
575	write(,) " inert_h2o_ice = ",inert_h2o_ice
576	! frost_albedo_threshold
577	write(,) "frost thickness threshold for albedo ?"
578	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
579	call getin("frost_albedo_threshold",
580	& frost_albedo_threshold)
581	write(,) " frost_albedo_threshold = ",
582	& frost_albedo_threshold
583
584	! call Titus crocus line -- DEFAULT IS NONE
585	write(,) "Titus crocus line ?"
586	tituscap=.false. ! default value
587	call getin("tituscap",tituscap)
588	write(,) "tituscap",tituscap
589
590
591	write(,) "photochemistry: include chemical species"
592	photochem=.false. ! default value
593	call getin("photochem",photochem)
594	write(,) " photochem = ",photochem
595
596
597	! THERMOSPHERE
598
599	write(,) "call thermosphere ?"
600	callthermos=.false. ! default value
601	call getin("callthermos",callthermos)
602	write(,) " callthermos = ",callthermos
603
604
605	write(,) " water included without cycle ",
606	& "(only if water=.false.)"
607	thermoswater=.false. ! default value
608	call getin("thermoswater",thermoswater)
609	write(,) " thermoswater = ",thermoswater
610
611	write(,) "call thermal conduction ?",
612	& " (only if callthermos=.true.)"
613	callconduct=.false. ! default value
614	call getin("callconduct",callconduct)
615	write(,) " callconduct = ",callconduct
616
617	write(,) "call EUV heating ?",
618	& " (only if callthermos=.true.)"
619	calleuv=.false. ! default value
620	call getin("calleuv",calleuv)
621	write(,) " calleuv = ",calleuv
622
623	write(,) "call molecular viscosity ?",
624	& " (only if callthermos=.true.)"
625	callmolvis=.false. ! default value
626	call getin("callmolvis",callmolvis)
627	write(,) " callmolvis = ",callmolvis
628
629	write(,) "call molecular diffusion ?",
630	& " (only if callthermos=.true.)"
631	callmoldiff=.false. ! default value
632	call getin("callmoldiff",callmoldiff)
633	write(,) " callmoldiff = ",callmoldiff
634
635
636	write(,) "call thermospheric photochemistry ?",
637	& " (only if callthermos=.true.)"
638	thermochem=.false. ! default value
639	call getin("thermochem",thermochem)
640	write(,) " thermochem = ",thermochem
641
642	write(,) "Method to include solar variability"
643	write(,) "0-> old method (using solarcondate); ",
644	& "1-> variability wit E10.7"
645	solvarmod=1
646	call getin("solvarmod",solvarmod)
647	write(,) " solvarmod = ",solvarmod
648
649	write(,) "date for solar flux calculation:",
650	& " (1985 < date < 2002)",
651	$ " (Only used if solvarmod=0)"
652	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
653	solarcondate=1993.4 ! default value
654	call getin("solarcondate",solarcondate)
655	write(,) " solarcondate = ",solarcondate
656
657	write(,) "Solar variability as observed for MY: "
658	write(,) "Only if solvarmod=1"
659	solvaryear=24
660	call getin("solvaryear",solvaryear)
661	write(,) " solvaryear = ",solvaryear
662
663	write(,) "UV heating efficiency:",
664	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
665	& "lower values may be used to compensate low 15 um cooling"
666	euveff=0.21 !default value
667	call getin("euveff",euveff)
668	write(,) " euveff = ", euveff
669
670	if (.not.callthermos) then
671	if (thermoswater) then
672	print*,'if thermoswater is set, callthermos must be true'
673	stop
674	endif
675	if (callconduct) then
676	print*,'if callconduct is set, callthermos must be true'
677	stop
678	endif
679	if (calleuv) then
680	print*,'if calleuv is set, callthermos must be true'
681	stop
682	endif
683	if (callmolvis) then
684	print*,'if callmolvis is set, callthermos must be true'
685	stop
686	endif
687	if (callmoldiff) then
688	print*,'if callmoldiff is set, callthermos must be true'
689	stop
690	endif
691	if (thermochem) then
692	print*,'if thermochem is set, callthermos must be true'
693	stop
694	endif
695	endif
696
697	! Test of incompatibility:
698	! if photochem is used, then water should be used too
699
700	if (photochem.and..not.water) then
701	print*,'if photochem is used, water should be used too'
702	stop
703	endif
704
705	! if callthermos is used, then thermoswater should be used too
706	! (if water not used already)
707
708	if (callthermos .and. .not.water) then
709	if (callthermos .and. .not.thermoswater) then
710	print*,'if callthermos is used, water or thermoswater
711	& should be used too'
712	stop
713	endif
714	endif
715
716	PRINT*,'--------------------------------------------'
717	PRINT*
718	PRINT*
719	ELSE
720	write(,)
721	write(,) 'Cannot read file callphys.def. Is it here ?'
722	stop
723	ENDIF
724
725	8000 FORMAT(t5,a12,l8)
726	8001 FORMAT(t5,a12,i8)
727
728	PRINT*
729	PRINT*,'inifis: daysec',daysec
730	PRINT*
731	PRINT*,'inifis: The radiative transfer is computed:'
732	PRINT*,' each ',iradia,' physical time-step'
733	PRINT,' or each ',iradiadtphys,' seconds'
734	PRINT*
735	! --------------------------------------------------------------
736	! Managing the Longwave radiative transfer
737	! --------------------------------------------------------------
738
739	! In most cases, the run just use the following values :
740	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
741	callemis=.true.
742	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
743	ilwd=1
744	ilwn=1 !2
745	ilwb=1 !2
746	linear=.true.
747	ncouche=3
748	alphan=0.4
749	semi=0
750
751	! BUT people working hard on the LW may want to read them in 'radia.def'
752	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
753
754	OPEN(99,file='radia.def',status='old',form='formatted'
755	. ,iostat=ierr)
756	IF(ierr.EQ.0) THEN
757	write(,) 'inifis: Reading radia.def !!!'
758	READ(99,fmt='(a)') ch1
759	READ(99,*) callemis
760	WRITE(*,8000) ch1,callemis
761
762	READ(99,fmt='(a)') ch1
763	READ(99,*) iradia
764	WRITE(*,8001) ch1,iradia
765
766	READ(99,fmt='(a)') ch1
767	READ(99,*) ilwd
768	WRITE(*,8001) ch1,ilwd
769
770	READ(99,fmt='(a)') ch1
771	READ(99,*) ilwn
772	WRITE(*,8001) ch1,ilwn
773
774	READ(99,fmt='(a)') ch1
775	READ(99,*) linear
776	WRITE(*,8000) ch1,linear
777
778	READ(99,fmt='(a)') ch1
779	READ(99,*) ncouche
780	WRITE(*,8001) ch1,ncouche
781
782	READ(99,fmt='(a)') ch1
783	READ(99,*) alphan
784	WRITE(,) ch1,alphan
785
786	READ(99,fmt='(a)') ch1
787	READ(99,*) ilwb
788	WRITE(*,8001) ch1,ilwb
789
790
791	READ(99,fmt='(a)') ch1
792	READ(99,'(l1)') callg2d
793	WRITE(*,8000) ch1,callg2d
794
795	READ(99,fmt='(a)') ch1
796	READ(99,*) semi
797	WRITE(,) ch1,semi
798	end if
799	CLOSE(99)
800
801	!-----------------------------------------------------------------------
802	! Some more initialization:
803	! ------------------------
804
805	! allocate "slope_mod" arrays
806	call ini_slope_mod(ngrid)
807
808	! allocate "comsaison_h" arrays
809	call ini_comsaison_h(ngrid)
810
811	! allocate "surfdat_h" arrays
812	call ini_surfdat_h(ngrid)
813
814	! allocate "comgeomfi_h" arrays
815	allocate(lati(ngrid))
816	allocate(long(ngrid))
817	allocate(area(ngrid))
818
819	! fill "comgeomfi_h" data
820	CALL SCOPY(ngrid,plon,1,long,1)
821	CALL SCOPY(ngrid,plat,1,lati,1)
822	CALL SCOPY(ngrid,parea,1,area,1)
823	totarea=SSUM(ngrid,area,1)
824
825	! allocate "comdiurn_h" data
826	allocate(sinlat(ngrid))
827	allocate(coslat(ngrid))
828	allocate(sinlon(ngrid))
829	allocate(coslon(ngrid))
830
831	! fill "comdiurn_h" data
832	DO ig=1,ngrid
833	sinlat(ig)=sin(plat(ig))
834	coslat(ig)=cos(plat(ig))
835	sinlon(ig)=sin(plon(ig))
836	coslon(ig)=cos(plon(ig))
837	ENDDO
838
839	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
840
841	! allocate "comsoil_h" arrays
842	call ini_comsoil_h(ngrid)
843
844	! set some variables in "dimradmars_mod"
845	call ini_dimradmars_mod(ngrid,nlayer)
846
847	! allocate arrays in "yomaer_h"
848	call ini_yomaer_h
849
850	! allocate arrays in "yomlw_h"
851	call ini_yomlw_h(ngrid)
852
853	! allocate arrays in "conc_mod"
854	call ini_conc_mod(ngrid,nlayer)
855
856	END

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

Download in other formats:

Original Format