Context Navigation

meso_inifis.F @ 226

Last change on this file since 226 was 226, checked in by aslmd, 13 years ago

MESOSCALE/LMDZ.MARS.new
--> modified to impact last changes

MESOSCALE/LMD_MM_MARS/makemeso
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq?.inc
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/call_meso_inifis?.inc
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
--> modified to get rid of ecri_phys

and make changes related to meso_physiq and meso_inifis

LMDZ.MARS/libf/phymars
--> see LMDZ.MARS/README

15/07/2011 == AS

Modified the mesoscale part so that the previous change by EM does not imply an error in the mesoscale case. More development is needed though to get the "varying z0" capability in the mesoscale model.
Worked on versions of meso_physiq and meso_inifis as close as possible to physiq and inifis for more continuity in the process of impacting changes (and even possibly to reach a common version of physiq and inifis).

The main point is to make the mesoscale significant specific parts

coded into include files in meso_inc so that meso_physiq and meso_inifis looks very close to physiq and inifis.

This is completely transparent for GCM users who does not need the

contents of meso_inc.
Slight cosmetic changes to physiq.f and inifis.F --- some of them e.g. to prepare convergence between meso_physiq and physiq

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 "meso_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	#ifndef MESOSCALE
141	datafile="/u/forget/WWW/datagcm/datafile"
142	#endif
143	call getin("datadir",datafile) ! default path
144	write(,) " datafile = ",trim(datafile)
145
146	write(,) "Run with or without tracer transport ?"
147	tracer=.false. ! default value
148	call getin("tracer",tracer)
149	write(,) " tracer = ",tracer
150
151	write(,) "Diurnal cycle ?"
152	write(,) "(if diurnal=False, diurnal averaged solar heating)"
153	diurnal=.true. ! default value
154	call getin("diurnal",diurnal)
155	write(,) " diurnal = ",diurnal
156
157	write(,) "Seasonal cycle ?"
158	write(,) "(if season=False, Ls stays constant, to value ",
159	& "set in 'start'"
160	season=.true. ! default value
161	call getin("season",season)
162	write(,) " season = ",season
163
164	write(,) "Write some extra output to the screen ?"
165	lwrite=.false. ! default value
166	call getin("lwrite",lwrite)
167	write(,) " lwrite = ",lwrite
168
169	write(,) "Save statistics in file stats.nc ?"
170	#ifdef MESOSCALE
171	callstats=.false. ! default value
172	#else
173	callstats=.true. ! default value
174	#endif
175	call getin("callstats",callstats)
176	write(,) " callstats = ",callstats
177
178	write(,) "Save EOF profiles in file 'profiles' for ",
179	& "Climate Database?"
180	calleofdump=.false. ! default value
181	call getin("calleofdump",calleofdump)
182	write(,) " calleofdump = ",calleofdump
183
184	write(,) "Dust scenario: 1=constant dust (read from startfi",
185	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
186	& "=4 Mars Year 24 from TES assimilation, ",
187	& "=24,25 or 26 :Mars Year 24,25 or 26 from TES assimilation"
188	iaervar=3 ! default value
189	call getin("iaervar",iaervar)
190	write(,) " iaervar = ",iaervar
191
192	write(,) "Reference (visible) dust opacity at 700 Pa ",
193	& "(matters only if iaervar=1)"
194	! NB: default value of tauvis is set/read in startfi.nc file
195	call getin("tauvis",tauvis)
196	write(,) " tauvis = ",tauvis
197
198	write(,) "Dust vertical distribution:"
199	write(,) "(=1 top set by topdustref parameter;",
200	& " =2 Viking scenario; =3 MGS scenario)"
201	iddist=3 ! default value
202	call getin("iddist",iddist)
203	write(,) " iddist = ",iddist
204
205	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
206	topdustref= 90.0 ! default value
207	call getin("topdustref",topdustref)
208	write(,) " topdustref = ",topdustref
209
210	write(,) "call radiative transfer ?"
211	callrad=.true. ! default value
212	call getin("callrad",callrad)
213	write(,) " callrad = ",callrad
214
215	write(,) "call NLTE radiative schemes ?",
216	& "(matters only if callrad=T)"
217	callnlte=.false. ! default value
218	call getin("callnlte",callnlte)
219	write(,) " callnlte = ",callnlte
220
221	write(,) "call CO2 NIR absorption ?",
222	& "(matters only if callrad=T)"
223	callnirco2=.false. ! default value
224	call getin("callnirco2",callnirco2)
225	write(,) " callnirco2 = ",callnirco2
226
227	write(,) "call turbulent vertical diffusion ?"
228	calldifv=.true. ! default value
229	call getin("calldifv",calldifv)
230	write(,) " calldifv = ",calldifv
231
232	write(,) "call thermals ?"
233	calltherm=.false. ! default value
234	call getin("calltherm",calltherm)
235	write(,) " calltherm = ",calltherm
236
237	write(,) "output thermal diagnostics ?"
238	outptherm=.false. ! default value
239	call getin("outptherm",outptherm)
240	write(,) " outptherm = ",outptherm
241
242	write(,) "call convective adjustment ?"
243	calladj=.true. ! default value
244	call getin("calladj",calladj)
245	write(,) " calladj = ",calladj
246
247	if (calltherm .and. (.not. calladj)) then
248	print*,'Convadj has to be activated when using thermals'
249	stop
250	endif
251
252	write(,) "call CO2 condensation ?"
253	callcond=.true. ! default value
254	call getin("callcond",callcond)
255	write(,) " callcond = ",callcond
256
257	write(,)"call thermal conduction in the soil ?"
258	callsoil=.true. ! default value
259	call getin("callsoil",callsoil)
260	write(,) " callsoil = ",callsoil
261
262
263	write(,)"call Lott's gravity wave/subgrid topography ",
264	& "scheme ?"
265	calllott=.true. ! default value
266	call getin("calllott",calllott)
267	write(,)" calllott = ",calllott
268
269
270	write(,)"rad.transfer is computed every iradia",
271	& " physical timestep"
272	iradia=1 ! default value
273	call getin("iradia",iradia)
274	write(,)" iradia = ",iradia
275
276
277	write(,)"Output of the exchange coefficient mattrix ?",
278	& "(for diagnostics only)"
279	callg2d=.false. ! default value
280	call getin("callg2d",callg2d)
281	write(,)" callg2d = ",callg2d
282
283	write(,)"Rayleigh scattering : (should be .false. for now)"
284	rayleigh=.false.
285	call getin("rayleigh",rayleigh)
286	write(,)" rayleigh = ",rayleigh
287
288
289	! TRACERS:
290
291	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
292	dustbin=0 ! default value
293	call getin("dustbin",dustbin)
294	write(,)" dustbin = ",dustbin
295
296	write(,)"Radiatively active dust ? (matters if dustbin>0)"
297	active=.false. ! default value
298	call getin("active",active)
299	write(,)" active = ",active
300
301	! Test of incompatibility:
302	! if active is used, then dustbin should be > 0
303
304	if (active.and.(dustbin.lt.1)) then
305	print*,'if active is used, then dustbin should > 0'
306	stop
307	endif
308
309	write(,)"use mass and number mixing ratios to predict",
310	& " dust size ?"
311	doubleq=.false. ! default value
312	call getin("doubleq",doubleq)
313	write(,)" doubleq = ",doubleq
314
315	submicron=.false. ! default value
316	call getin("submicron",submicron)
317	write(,)" submicron = ",submicron
318
319	! Test of incompatibility:
320	! if doubleq is used, then dustbin should be 2
321
322	if (doubleq.and.(dustbin.ne.2)) then
323	print*,'if doubleq is used, then dustbin should be 2'
324	stop
325	endif
326	if (doubleq.and.submicron.and.(nqmx.LT.3)) then
327	print*,'If doubleq is used with a submicron tracer,'
328	print*,' then the number of tracers has to be'
329	print*,' larger than 3.'
330	stop
331	endif
332
333	write(,)"dust lifted by GCM surface winds ?"
334	lifting=.false. ! default value
335	call getin("lifting",lifting)
336	write(,)" lifting = ",lifting
337
338	! Test of incompatibility:
339	! if lifting is used, then dustbin should be > 0
340
341	if (lifting.and.(dustbin.lt.1)) then
342	print*,'if lifting is used, then dustbin should > 0'
343	stop
344	endif
345
346	write(,)" dust lifted by dust devils ?"
347	callddevil=.false. !default value
348	call getin("callddevil",callddevil)
349	write(,)" callddevil = ",callddevil
350
351
352	! Test of incompatibility:
353	! if dustdevil is used, then dustbin should be > 0
354
355	if (callddevil.and.(dustbin.lt.1)) then
356	print*,'if dustdevil is used, then dustbin should > 0'
357	stop
358	endif
359
360	write(,)"Dust scavenging by CO2 snowfall ?"
361	scavenging=.false. ! default value
362	call getin("scavenging",scavenging)
363	write(,)" scavenging = ",scavenging
364
365
366	! Test of incompatibility:
367	! if scavenging is used, then dustbin should be > 0
368
369	if (scavenging.and.(dustbin.lt.1)) then
370	print*,'if scavenging is used, then dustbin should > 0'
371	stop
372	endif
373
374	write(,) "Gravitationnal sedimentation ?"
375	sedimentation=.true. ! default value
376	call getin("sedimentation",sedimentation)
377	write(,) " sedimentation = ",sedimentation
378
379	write(,) "Radiatively active transported atmospheric ",
380	& "water ice ?"
381	activice=.false. ! default value
382	call getin("activice",activice)
383	write(,) " activice = ",activice
384
385	write(,) "Compute water cycle ?"
386	water=.false. ! default value
387	call getin("water",water)
388	write(,) " water = ",water
389
390	! Test of incompatibility:
391
392	if (activice.and..not.water) then
393	print*,'if activice is used, water should be used too'
394	stop
395	endif
396
397	if (water.and..not.tracer) then
398	print*,'if water is used, tracer should be used too'
399	stop
400	endif
401
402	! Test of incompatibility:
403
404	write(,) "Permanent water caps at poles ?",
405	& " .true. is RECOMMENDED"
406	write(,) "(with .true., North cap is a source of water ",
407	& "and South pole is a cold trap)"
408	caps=.true. ! default value
409	call getin("caps",caps)
410	write(,) " caps = ",caps
411
412	write(,) "photochemistry: include chemical species"
413	photochem=.false. ! default value
414	call getin("photochem",photochem)
415	write(,) " photochem = ",photochem
416
417
418	! THERMOSPHERE
419
420	write(,) "call thermosphere ?"
421	callthermos=.false. ! default value
422	call getin("callthermos",callthermos)
423	write(,) " callthermos = ",callthermos
424
425
426	write(,) " water included without cycle ",
427	& "(only if water=.false.)"
428	thermoswater=.false. ! default value
429	call getin("thermoswater",thermoswater)
430	write(,) " thermoswater = ",thermoswater
431
432	write(,) "call thermal conduction ?",
433	& " (only if callthermos=.true.)"
434	callconduct=.false. ! default value
435	call getin("callconduct",callconduct)
436	write(,) " callconduct = ",callconduct
437
438	write(,) "call EUV heating ?",
439	& " (only if callthermos=.true.)"
440	calleuv=.false. ! default value
441	call getin("calleuv",calleuv)
442	write(,) " calleuv = ",calleuv
443
444	write(,) "call molecular viscosity ?",
445	& " (only if callthermos=.true.)"
446	callmolvis=.false. ! default value
447	call getin("callmolvis",callmolvis)
448	write(,) " callmolvis = ",callmolvis
449
450	write(,) "call molecular diffusion ?",
451	& " (only if callthermos=.true.)"
452	callmoldiff=.false. ! default value
453	call getin("callmoldiff",callmoldiff)
454	write(,) " callmoldiff = ",callmoldiff
455
456
457	write(,) "call thermospheric photochemistry ?",
458	& " (only if callthermos=.true.)"
459	thermochem=.false. ! default value
460	call getin("thermochem",thermochem)
461	write(,) " thermochem = ",thermochem
462
463	write(,) "date for solar flux calculation:",
464	& " (1985 < date < 2002)"
465	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
466	solarcondate=1993.4 ! default value
467	call getin("solarcondate",solarcondate)
468	write(,) " solarcondate = ",solarcondate
469
470
471	if (.not.callthermos) then
472	if (thermoswater) then
473	print*,'if thermoswater is set, callthermos must be true'
474	stop
475	endif
476	if (callconduct) then
477	print*,'if callconduct is set, callthermos must be true'
478	stop
479	endif
480	if (calleuv) then
481	print*,'if calleuv is set, callthermos must be true'
482	stop
483	endif
484	if (callmolvis) then
485	print*,'if callmolvis is set, callthermos must be true'
486	stop
487	endif
488	if (callmoldiff) then
489	print*,'if callmoldiff is set, callthermos must be true'
490	stop
491	endif
492	if (thermochem) then
493	print*,'if thermochem is set, callthermos must be true'
494	stop
495	endif
496	endif
497
498	! Test of incompatibility:
499	! if photochem is used, then water should be used too
500
501	if (photochem.and..not.water) then
502	print*,'if photochem is used, water should be used too'
503	stop
504	endif
505
506	! if callthermos is used, then thermoswater should be used too
507	! (if water not used already)
508
509	if (callthermos .and. .not.water) then
510	if (callthermos .and. .not.thermoswater) then
511	print*,'if callthermos is used, water or thermoswater
512	& should be used too'
513	stop
514	endif
515	endif
516
517	PRINT*,'--------------------------------------------'
518	PRINT*
519	PRINT*
520	ELSE
521	write(,)
522	write(,) 'Cannot read file callphys.def. Is it here ?'
523	stop
524	ENDIF
525
526	8000 FORMAT(t5,a12,l8)
527	8001 FORMAT(t5,a12,i8)
528
529	PRINT*
530	PRINT*,'inifis: daysec',daysec
531	PRINT*
532	PRINT*,'inifis: The radiative transfer is computed:'
533	PRINT*,' each ',iradia,' physical time-step'
534	PRINT,' or each ',iradiadtphys,' seconds'
535	PRINT*
536	! --------------------------------------------------------------
537	! Managing the Longwave radiative transfer
538	! --------------------------------------------------------------
539
540	! In most cases, the run just use the following values :
541	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
542	callemis=.true.
543	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
544	ilwd=1
545	ilwn=1 !2
546	ilwb=1 !2
547	linear=.true.
548	ncouche=3
549	alphan=0.4
550	semi=0
551
552	! BUT people working hard on the LW may want to read them in 'radia.def'
553	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
554
555	OPEN(99,file='radia.def',status='old',form='formatted'
556	. ,iostat=ierr)
557	IF(ierr.EQ.0) THEN
558	write(,) 'inifis: Reading radia.def !!!'
559	READ(99,fmt='(a)') ch1
560	READ(99,*) callemis
561	WRITE(*,8000) ch1,callemis
562
563	READ(99,fmt='(a)') ch1
564	READ(99,*) iradia
565	WRITE(*,8001) ch1,iradia
566
567	READ(99,fmt='(a)') ch1
568	READ(99,*) ilwd
569	WRITE(*,8001) ch1,ilwd
570
571	READ(99,fmt='(a)') ch1
572	READ(99,*) ilwn
573	WRITE(*,8001) ch1,ilwn
574
575	READ(99,fmt='(a)') ch1
576	READ(99,*) linear
577	WRITE(*,8000) ch1,linear
578
579	READ(99,fmt='(a)') ch1
580	READ(99,*) ncouche
581	WRITE(*,8001) ch1,ncouche
582
583	READ(99,fmt='(a)') ch1
584	READ(99,*) alphan
585	WRITE(,) ch1,alphan
586
587	READ(99,fmt='(a)') ch1
588	READ(99,*) ilwb
589	WRITE(*,8001) ch1,ilwb
590
591
592	READ(99,fmt='(a)') ch1
593	READ(99,'(l1)') callg2d
594	WRITE(*,8000) ch1,callg2d
595
596	READ(99,fmt='(a)') ch1
597	READ(99,*) semi
598	WRITE(,) ch1,semi
599	end if
600	CLOSE(99)
601
602	!-----------------------------------------------------------------------
603	! Some more initialization:
604	! ------------------------
605
606	! in 'comgeomfi.h'
607	CALL SCOPY(ngrid,plon,1,long,1)
608	CALL SCOPY(ngrid,plat,1,lati,1)
609	CALL SCOPY(ngrid,parea,1,area,1)
610	totarea=SSUM(ngridmx,area,1)
611
612	! in 'comdiurn.h'
613	DO ig=1,ngrid
614	sinlat(ig)=sin(plat(ig))
615	coslat(ig)=cos(plat(ig))
616	sinlon(ig)=sin(plon(ig))
617	coslon(ig)=cos(plon(ig))
618	ENDDO
619
620	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
621
622	! managing the tracers, and tests:
623	! -------------------------------
624	! Ehouarn: removed; as these tests are now done in initracer.F
625	! if(tracer) then
626	!
627	!! when photochem is used, nqchem_min is the rank
628	!! of the first chemical species
629	!
630	!! Ehouarn: nqchem_min is now meaningless and no longer used
631	!! nqchem_min = 1
632	! if (photochem .or. callthermos) then
633	! chem = .true.
634	! end if
635	!
636	! if (water .or. thermoswater) h2o = .true.
637	!
638	!! TESTS
639	!
640	! print*,'inifis: TRACERS:'
641	! write(,) " chem=",chem," h2o=",h2o
642	!! write(,) " doubleq=",doubleq
643	!! write(,) " dustbin=",dustbin
644	!
645	! if ((doubleq).and.(h2o).and.
646	! $ (chem)) then
647	! print*,' 2 dust tracers (doubleq)'
648	! print*,' 1 water vapour tracer'
649	! print*,' 1 water ice tracer'
650	! print*,nqmx-4,' chemistry tracers'
651	! endif
652	!
653	! if ((doubleq).and.(h2o).and.
654	! $ .not.(chem)) then
655	! print*,' 2 dust tracers (doubleq)'
656	! print*,' 1 water vapour tracer'
657	! print*,' 1 water ice tracer'
658	! if (nqmx.LT.4) then
659	! print*,'nqmx should be at least equal to'
660	! print*,'4 with these options.'
661	! stop
662	! endif
663	! endif
664	!
665	! if (.not.(doubleq).and.(h2o).and.
666	! $ (chem)) then
667	! if (dustbin.gt.0) then
668	! print*,dustbin,' dust bins'
669	! endif
670	! print*,nqmx-2-dustbin,' chemistry tracers'
671	! print*,' 1 water vapour tracer'
672	! print*,' 1 water ice tracer'
673	! endif
674	!
675	! if (.not.(doubleq).and.(h2o).and.
676	! $ .not.(chem)) then
677	! if (dustbin.gt.0) then
678	! print*,dustbin,' dust bins'
679	! endif
680	! print*,' 1 water vapour tracer'
681	! print*,' 1 water ice tracer'
682	! if (nqmx.gt.(dustbin+2)) then
683	! print*,'nqmx should be ',(dustbin+2),
684	! $ ' with these options...'
685	! print*,'(or check callphys.def)'
686	! endif
687	! if (nqmx.lt.(dustbin+2)) then
688	! write(,) "inifis: nqmx.lt.(dustbin+2)"
689	! stop
690	! endif
691	! endif
692	!
693	! endif ! of if (tracer)
694	!
695	! RETURN
696	END

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

Download in other formats:

Original Format