Context Navigation

inifis.F @ 146

Last change on this file since 146 was 79, checked in by emillour, 14 years ago
Mars GCM: small correction in inifis.F and datafile.h EM
File size: 20.3 KB

Line
1	SUBROUTINE inifis(ngrid,nlayer,
2	$ day_ini,pdaysec,ptimestep,
3	$ plat,plon,parea,
4	$ prad,pg,pr,pcpp)
5	!
6	!=======================================================================
7	!
8	! purpose:
9	! -------
10	!
11	! Initialisation for the physical parametrisations of the LMD
12	! martian atmospheric general circulation modele.
13	!
14	! author: Frederic Hourdin 15 / 10 /93
15	! -------
16	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
17	! Ehouarn Millour (oct. 2008) tracers are now identified
18	! by their names and may not be contiguously
19	! stored in the q(:,:,:,:) array
20	! E.M. (june 2009) use getin routine to load parameters
21	!
22	!
23	! arguments:
24	! ----------
25	!
26	! input:
27	! ------
28	!
29	! ngrid Size of the horizontal grid.
30	! All internal loops are performed on that grid.
31	! nlayer Number of vertical layers.
32	! pdayref Day of reference for the simulation
33	! pday Number of days counted from the North. Spring
34	! equinoxe.
35	!
36	!=======================================================================
37	!
38	!-----------------------------------------------------------------------
39	! declarations:
40	! -------------
41	! to use 'getin'
42	USE ioipsl_getincom
43	IMPLICIT NONE
44	#include "dimensions.h"
45	#include "dimphys.h"
46	#include "planete.h"
47	#include "comcstfi.h"
48	#include "comsaison.h"
49	#include "comdiurn.h"
50	#include "comgeomfi.h"
51	#include "callkeys.h"
52	#include "surfdat.h"
53	#include "dimradmars.h"
54	#include "yomaer.h"
55	#include "datafile.h"
56
57	REAL prad,pg,pr,pcpp,pdaysec,ptimestep
58
59	INTEGER ngrid,nlayer
60	REAL plat(ngrid),plon(ngrid),parea(ngridmx)
61	integer day_ini
62	INTEGER ig,ierr
63
64	! EXTERNAL iniorbit,orbite
65	EXTERNAL SSUM
66	REAL SSUM
67
68	CHARACTER ch1*12
69	CHARACTER ch80*80
70
71	! logical chem, h2o
72
73	! chem = .false.
74	! h2o = .false.
75
76	rad=prad
77	daysec=pdaysec
78	dtphys=ptimestep
79	cpp=pcpp
80	g=pg
81	r=pr
82	rcp=r/cpp
83
84	! --------------------------------------------------------
85	! The usual Tests
86	! --------------
87	IF (nlayer.NE.nlayermx) THEN
88	PRINT*,'STOP in inifis'
89	PRINT*,'Probleme de dimensions :'
90	PRINT*,'nlayer = ',nlayer
91	PRINT*,'nlayermx = ',nlayermx
92	STOP
93	ENDIF
94
95	IF (ngrid.NE.ngridmx) THEN
96	PRINT*,'STOP in inifis'
97	PRINT*,'Probleme de dimensions :'
98	PRINT*,'ngrid = ',ngrid
99	PRINT*,'ngridmx = ',ngridmx
100	STOP
101	ENDIF
102
103	! --------------------------------------------------------------
104	! Reading the "callphys.def" file controlling some key options
105	! --------------------------------------------------------------
106
107	! check that 'callphys.def' file is around
108	OPEN(99,file='callphys.def',status='old',form='formatted'
109	& ,iostat=ierr)
110	CLOSE(99)
111
112	IF(ierr.EQ.0) THEN
113	PRINT*
114	PRINT*
115	PRINT*,'--------------------------------------------'
116	PRINT*,' inifis: Parameters for the physics (callphys.def)'
117	PRINT*,'--------------------------------------------'
118
119	write(,) "Directory where external input files are:"
120	datafile="/u/forget/WWW/datagcm/datafile"
121	call getin("datadir",datafile) ! default path
122	write(,) " datafile = ",trim(datafile)
123
124	write(,) "Run with or without tracer transport ?"
125	tracer=.false. ! default value
126	call getin("tracer",tracer)
127	write(,) " tracer = ",tracer
128
129	write(,) "Diurnal cycle ?"
130	write(,) "(if diurnal=False, diurnal averaged solar heating)"
131	diurnal=.true. ! default value
132	call getin("diurnal",diurnal)
133	write(,) " diurnal = ",diurnal
134
135	write(,) "Seasonal cycle ?"
136	write(,) "(if season=False, Ls stays constant, to value ",
137	& "set in 'start'"
138	season=.true. ! default value
139	call getin("season",season)
140	write(,) " season = ",season
141
142	write(,) "Write some extra output to the screen ?"
143	lwrite=.false. ! default value
144	call getin("lwrite",lwrite)
145	write(,) " lwrite = ",lwrite
146
147	write(,) "Save statistics in file stats.nc ?"
148	callstats=.true. ! default value
149	call getin("callstats",callstats)
150	write(,) " callstats = ",callstats
151
152	write(,) "Save EOF profiles in file 'profiles' for ",
153	& "Climate Database?"
154	calleofdump=.false. ! default value
155	call getin("calleofdump",calleofdump)
156	write(,) " calleofdump = ",calleofdump
157
158	write(,) "Dust scenario: 1=constant dust (read from startfi",
159	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
160	& "=4 Mars Year 24 from TES assimilation, ",
161	& "=24,25 or 26 :Mars Year 24,25 or 26 from TES assimilation"
162	iaervar=3 ! default value
163	call getin("iaervar",iaervar)
164	write(,) " iaervar = ",iaervar
165
166	write(,) "Reference (visible) dust opacity at 700 Pa ",
167	& "(matters only if iaervar=1)"
168	! NB: default value of tauvis is set/read in startfi.nc file
169	call getin("tauvis",tauvis)
170	write(,) " tauvis = ",tauvis
171
172	write(,) "Dust vertical distribution:"
173	write(,) "(=1 top set by topdustref parameter;",
174	& " =2 Viking scenario; =3 MGS scenario)"
175	iddist=3 ! default value
176	call getin("iddist",iddist)
177	write(,) " iddist = ",iddist
178
179	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
180	topdustref= 90.0 ! default value
181	call getin("topdustref",topdustref)
182	write(,) " topdustref = ",topdustref
183
184	write(,) "call radiative transfer ?"
185	callrad=.true. ! default value
186	call getin("callrad",callrad)
187	write(,) " callrad = ",callrad
188
189	write(,) "call NLTE radiative schemes ?",
190	& "(matters only if callrad=T)"
191	callnlte=.false. ! default value
192	call getin("callnlte",callnlte)
193	write(,) " callnlte = ",callnlte
194
195	write(,) "call CO2 NIR absorption ?",
196	& "(matters only if callrad=T)"
197	callnirco2=.false. ! default value
198	call getin("callnirco2",callnirco2)
199	write(,) " callnirco2 = ",callnirco2
200
201	write(,) "call turbulent vertical diffusion ?"
202	calldifv=.true. ! default value
203	call getin("calldifv",calldifv)
204	write(,) " calldifv = ",calldifv
205
206	write(,) "call convective adjustment ?"
207	calladj=.true. ! default value
208	call getin("calladj",calladj)
209	write(,) " calladj = ",calladj
210
211
212	write(,) "call CO2 condensation ?"
213	callcond=.true. ! default value
214	call getin("callcond",callcond)
215	write(,) " callcond = ",callcond
216
217	write(,)"call thermal conduction in the soil ?"
218	callsoil=.true. ! default value
219	call getin("callsoil",callsoil)
220	write(,) " callsoil = ",callsoil
221
222
223	write(,)"call Lott's gravity wave/subgrid topography ",
224	& "scheme ?"
225	calllott=.true. ! default value
226	call getin("calllott",calllott)
227	write(,)" calllott = ",calllott
228
229
230	write(,)"rad.transfer is computed every iradia",
231	& " physical timestep"
232	iradia=1 ! default value
233	call getin("iradia",iradia)
234	write(,)" iradia = ",iradia
235
236
237	write(,)"Output of the exchange coefficient mattrix ?",
238	& "(for diagnostics only)"
239	callg2d=.false. ! default value
240	call getin("callg2d",callg2d)
241	write(,)" callg2d = ",callg2d
242
243	write(,)"Rayleigh scattering : (should be .false. for now)"
244	rayleigh=.false.
245	call getin("rayleigh",rayleigh)
246	write(,)" rayleigh = ",rayleigh
247
248
249	! TRACERS:
250
251	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
252	dustbin=0 ! default value
253	call getin("dustbin",dustbin)
254	write(,)" dustbin = ",dustbin
255
256	write(,)"Radiatively active dust ? (matters if dustbin>0)"
257	active=.false. ! default value
258	call getin("active",active)
259	write(,)" active = ",active
260
261	! Test of incompatibility:
262	! if active is used, then dustbin should be > 0
263
264	if (active.and.(dustbin.lt.1)) then
265	print*,'if active is used, then dustbin should > 0'
266	stop
267	endif
268
269	write(,)"use mass and number mixing ratios to predict",
270	& " dust size ?"
271	doubleq=.false. ! default value
272	call getin("doubleq",doubleq)
273	write(,)" doubleq = ",doubleq
274
275	submicron=.false. ! default value
276	call getin("submicron",submicron)
277	write(,)" submicron = ",submicron
278
279	! Test of incompatibility:
280	! if doubleq is used, then dustbin should be 2
281
282	if (doubleq.and.(dustbin.ne.2)) then
283	print*,'if doubleq is used, then dustbin should be 2'
284	stop
285	endif
286	if (doubleq.and.submicron.and.(nqmx.LT.3)) then
287	print*,'If doubleq is used with a submicron tracer,'
288	print*,' then the number of tracers has to be'
289	print*,' larger than 3.'
290	stop
291	endif
292
293	write(,)"dust lifted by GCM surface winds ?"
294	lifting=.false. ! default value
295	call getin("lifting",lifting)
296	write(,)" lifting = ",lifting
297
298	! Test of incompatibility:
299	! if lifting is used, then dustbin should be > 0
300
301	if (lifting.and.(dustbin.lt.1)) then
302	print*,'if lifting is used, then dustbin should > 0'
303	stop
304	endif
305
306	write(,)" dust lifted by dust devils ?"
307	callddevil=.false. !default value
308	call getin("callddevil",callddevil)
309	write(,)" callddevil = ",callddevil
310
311
312	! Test of incompatibility:
313	! if dustdevil is used, then dustbin should be > 0
314
315	if (callddevil.and.(dustbin.lt.1)) then
316	print*,'if dustdevil is used, then dustbin should > 0'
317	stop
318	endif
319
320	write(,)"Dust scavenging by CO2 snowfall ?"
321	scavenging=.false. ! default value
322	call getin("scavenging",scavenging)
323	write(,)" scavenging = ",scavenging
324
325
326	! Test of incompatibility:
327	! if scavenging is used, then dustbin should be > 0
328
329	if (scavenging.and.(dustbin.lt.1)) then
330	print*,'if scavenging is used, then dustbin should > 0'
331	stop
332	endif
333
334	write(,) "Gravitationnal sedimentation ?"
335	sedimentation=.true. ! default value
336	call getin("sedimentation",sedimentation)
337	write(,) " sedimentation = ",sedimentation
338
339	write(,) "Radiatively active transported atmospheric ",
340	& "water ice ?"
341	activice=.false. ! default value
342	call getin("activice",activice)
343	write(,) " activice = ",activice
344
345	write(,) "Compute water cycle ?"
346	water=.false. ! default value
347	call getin("water",water)
348	write(,) " water = ",water
349
350	! Test of incompatibility:
351
352	if (activice.and..not.water) then
353	print*,'if activice is used, water should be used too'
354	stop
355	endif
356
357	if (water.and..not.tracer) then
358	print*,'if water is used, tracer should be used too'
359	stop
360	endif
361
362	! Test of incompatibility:
363
364	write(,) "Permanent water caps at poles ?",
365	& " .true. is RECOMMENDED"
366	write(,) "(with .true., North cap is a source of water ",
367	& "and South pole is a cold trap)"
368	caps=.true. ! default value
369	call getin("caps",caps)
370	write(,) " caps = ",caps
371
372	write(,) "photochemistry: include chemical species"
373	photochem=.false. ! default value
374	call getin("photochem",photochem)
375	write(,) " photochem = ",photochem
376
377
378	! THERMOSPHERE
379
380	write(,) "call thermosphere ?"
381	callthermos=.false. ! default value
382	call getin("callthermos",callthermos)
383	write(,) " callthermos = ",callthermos
384
385
386	write(,) " water included without cycle ",
387	& "(only if water=.false.)"
388	thermoswater=.false. ! default value
389	call getin("thermoswater",thermoswater)
390	write(,) " thermoswater = ",thermoswater
391
392	write(,) "call thermal conduction ?",
393	& " (only if callthermos=.true.)"
394	callconduct=.false. ! default value
395	call getin("callconduct",callconduct)
396	write(,) " callconduct = ",callconduct
397
398	write(,) "call EUV heating ?",
399	& " (only if callthermos=.true.)"
400	calleuv=.false. ! default value
401	call getin("calleuv",calleuv)
402	write(,) " calleuv = ",calleuv
403
404	write(,) "call molecular viscosity ?",
405	& " (only if callthermos=.true.)"
406	callmolvis=.false. ! default value
407	call getin("callmolvis",callmolvis)
408	write(,) " callmolvis = ",callmolvis
409
410	write(,) "call molecular diffusion ?",
411	& " (only if callthermos=.true.)"
412	callmoldiff=.false. ! default value
413	call getin("callmoldiff",callmoldiff)
414	write(,) " callmoldiff = ",callmoldiff
415
416
417	write(,) "call thermospheric photochemistry ?",
418	& " (only if callthermos=.true.)"
419	thermochem=.false. ! default value
420	call getin("thermochem",thermochem)
421	write(,) " thermochem = ",thermochem
422
423	write(,) "date for solar flux calculation:",
424	& " (1985 < date < 2002)"
425	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
426	solarcondate=1993.4 ! default value
427	call getin("solarcondate",solarcondate)
428	write(,) " solarcondate = ",solarcondate
429
430
431	if (.not.callthermos) then
432	if (thermoswater) then
433	print*,'if thermoswater is set, callthermos must be true'
434	stop
435	endif
436	if (callconduct) then
437	print*,'if callconduct is set, callthermos must be true'
438	stop
439	endif
440	if (calleuv) then
441	print*,'if calleuv is set, callthermos must be true'
442	stop
443	endif
444	if (callmolvis) then
445	print*,'if callmolvis is set, callthermos must be true'
446	stop
447	endif
448	if (callmoldiff) then
449	print*,'if callmoldiff is set, callthermos must be true'
450	stop
451	endif
452	if (thermochem) then
453	print*,'if thermochem is set, callthermos must be true'
454	stop
455	endif
456	endif
457
458	! Test of incompatibility:
459	! if photochem is used, then water should be used too
460
461	if (photochem.and..not.water) then
462	print*,'if photochem is used, water should be used too'
463	stop
464	endif
465
466	! if callthermos is used, then thermoswater should be used too
467	! (if water not used already)
468
469	if (callthermos .and. .not.water) then
470	if (callthermos .and. .not.thermoswater) then
471	print*,'if callthermos is used, water or thermoswater
472	& should be used too'
473	stop
474	endif
475	endif
476
477	PRINT*,'--------------------------------------------'
478	PRINT*
479	PRINT*
480	ELSE
481	write(,)
482	write(,) 'Cannot read file callphys.def. Is it here ?'
483	stop
484	ENDIF
485
486	8000 FORMAT(t5,a12,l8)
487	8001 FORMAT(t5,a12,i8)
488
489	PRINT*
490	PRINT*,'inifis: daysec',daysec
491	PRINT*
492	PRINT*,'inifis: The radiative transfer is computed:'
493	PRINT*,' each ',iradia,' physical time-step'
494	PRINT,' or each ',iradiadtphys,' seconds'
495	PRINT*
496	! --------------------------------------------------------------
497	! Managing the Longwave radiative transfer
498	! --------------------------------------------------------------
499
500	! In most cases, the run just use the following values :
501	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
502	callemis=.true.
503	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
504	ilwd=1
505	ilwn=1 !2
506	ilwb=1 !2
507	linear=.true.
508	ncouche=3
509	alphan=0.4
510	semi=0
511
512	! BUT people working hard on the LW may want to read them in 'radia.def'
513	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
514
515	OPEN(99,file='radia.def',status='old',form='formatted'
516	. ,iostat=ierr)
517	IF(ierr.EQ.0) THEN
518	write(,) 'inifis: Reading radia.def !!!'
519	READ(99,fmt='(a)') ch1
520	READ(99,*) callemis
521	WRITE(*,8000) ch1,callemis
522
523	READ(99,fmt='(a)') ch1
524	READ(99,*) iradia
525	WRITE(*,8001) ch1,iradia
526
527	READ(99,fmt='(a)') ch1
528	READ(99,*) ilwd
529	WRITE(*,8001) ch1,ilwd
530
531	READ(99,fmt='(a)') ch1
532	READ(99,*) ilwn
533	WRITE(*,8001) ch1,ilwn
534
535	READ(99,fmt='(a)') ch1
536	READ(99,*) linear
537	WRITE(*,8000) ch1,linear
538
539	READ(99,fmt='(a)') ch1
540	READ(99,*) ncouche
541	WRITE(*,8001) ch1,ncouche
542
543	READ(99,fmt='(a)') ch1
544	READ(99,*) alphan
545	WRITE(,) ch1,alphan
546
547	READ(99,fmt='(a)') ch1
548	READ(99,*) ilwb
549	WRITE(*,8001) ch1,ilwb
550
551
552	READ(99,fmt='(a)') ch1
553	READ(99,'(l1)') callg2d
554	WRITE(*,8000) ch1,callg2d
555
556	READ(99,fmt='(a)') ch1
557	READ(99,*) semi
558	WRITE(,) ch1,semi
559	end if
560	CLOSE(99)
561
562	!-----------------------------------------------------------------------
563	! Some more initialization:
564	! ------------------------
565
566	CALL SCOPY(ngrid,plon,1,long,1)
567	CALL SCOPY(ngrid,plat,1,lati,1)
568	CALL SCOPY(ngrid,parea,1,area,1)
569	totarea=SSUM(ngridmx,area,1)
570
571	DO ig=1,ngrid
572	sinlat(ig)=sin(plat(ig))
573	coslat(ig)=cos(plat(ig))
574	sinlon(ig)=sin(plon(ig))
575	coslon(ig)=cos(plon(ig))
576	ENDDO
577
578	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
579
580	! managing the tracers, and tests:
581	! -------------------------------
582	! Ehouarn: removed; as these tests are now done in initracer.F
583	! if(tracer) then
584	!
585	!! when photochem is used, nqchem_min is the rank
586	!! of the first chemical species
587	!
588	!! Ehouarn: nqchem_min is now meaningless and no longer used
589	!! nqchem_min = 1
590	! if (photochem .or. callthermos) then
591	! chem = .true.
592	! end if
593	!
594	! if (water .or. thermoswater) h2o = .true.
595	!
596	!! TESTS
597	!
598	! print*,'inifis: TRACERS:'
599	! write(,) " chem=",chem," h2o=",h2o
600	!! write(,) " doubleq=",doubleq
601	!! write(,) " dustbin=",dustbin
602	!
603	! if ((doubleq).and.(h2o).and.
604	! $ (chem)) then
605	! print*,' 2 dust tracers (doubleq)'
606	! print*,' 1 water vapour tracer'
607	! print*,' 1 water ice tracer'
608	! print*,nqmx-4,' chemistry tracers'
609	! endif
610	!
611	! if ((doubleq).and.(h2o).and.
612	! $ .not.(chem)) then
613	! print*,' 2 dust tracers (doubleq)'
614	! print*,' 1 water vapour tracer'
615	! print*,' 1 water ice tracer'
616	! if (nqmx.LT.4) then
617	! print*,'nqmx should be at least equal to'
618	! print*,'4 with these options.'
619	! stop
620	! endif
621	! endif
622	!
623	! if (.not.(doubleq).and.(h2o).and.
624	! $ (chem)) then
625	! if (dustbin.gt.0) then
626	! print*,dustbin,' dust bins'
627	! endif
628	! print*,nqmx-2-dustbin,' chemistry tracers'
629	! print*,' 1 water vapour tracer'
630	! print*,' 1 water ice tracer'
631	! endif
632	!
633	! if (.not.(doubleq).and.(h2o).and.
634	! $ .not.(chem)) then
635	! if (dustbin.gt.0) then
636	! print*,dustbin,' dust bins'
637	! endif
638	! print*,' 1 water vapour tracer'
639	! print*,' 1 water ice tracer'
640	! if (nqmx.gt.(dustbin+2)) then
641	! print*,'nqmx should be ',(dustbin+2),
642	! $ ' with these options...'
643	! print*,'(or check callphys.def)'
644	! endif
645	! if (nqmx.lt.(dustbin+2)) then
646	! write(,) "inifis: nqmx.lt.(dustbin+2)"
647	! stop
648	! endif
649	! endif
650	!
651	! endif ! of if (tracer)
652	!
653	! RETURN
654	END

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

Download in other formats:

Original Format