Context Navigation

meso_inifis.F @ 228

Last change on this file since 228 was 228, checked in by aslmd, 14 years ago
MESOSCALE: corrected the diagfi bug. updated test case and check those are still working OK. other changes are cosmetic
File size: 21.4 KB

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

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

Download in other formats:

Original Format