Context Navigation

inifis.F @ 161

Last change on this file since 161 was 161, checked in by acolaitis, 13 years ago

================================================
======== IMPLEMENTATION OF THERMALS ============
================================================

Author: A. Colaitis (2011-06-16)

The main goal of this revision is to start including the thermals into the model
for development purposes. Users should not use the thermals yet, as
several major configuration changes still need to be done.

This version includes :

updraft and downdraft parametrizations
velocity in the thermal, including drag
plume height analysis
closure equation
updraft transport of heat, tracers and momentum
downdraft transport of heat

This model should not be used without upcoming developments, namely :

downdraft transport of tracers and momentum
updraft & downdraft transport of q2 (tke)
revision of vdif_kc to compute q2 for non-stratified cases

Thermals could also include in a later revision :

momentum loss during transport (horizontal drag)

Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90

================================================
================================================

M libf/phymars/callkeys.h
M libf/phymars/inifis.F

Added new control flags to call the thermals :

calltherm (false by default) <- to call thermals
outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)

================================================

M libf/phymars/vdifc.F
^{------> added a temporary output for thermal-related diagnostics}

M libf/phymars/testphys1d.F
^{------> added treatment for a initialization from a profile of neutral gas (ar)}

-> will be transformed in a decaying tracer for thermal diagnostics

M libf/phymars/physiq.F
^{------> added a section to call the thermals}

-> changed the call to convadj
-> added thermal-related outputs for diagnostics

M libf/phymars/convadj.F
^{------> takes now into account the height of thermals to execute convective adjustment}

=> note : convective adjustment needs to be activated when using thermals, in case of a

second instable layer above the thermals

================================================

A libf/phymars/calltherm_interface.F90
^{------> Interface between physiq.F and the thermals}

A libf/phymars/calltherm_mars.F90
^{------> Routine running the sub-timestep of the thermals}

A libf/phymars/thermcell_main_mars.F90
^{------> Main thermals routine specific to Martian physics}

A libf/phymars/thermcell_dqupdown.F90
^{------> Thermals subroutine computing transport of quantities by updrafts and downdrafts}

A libf/phymars/thermcell.F90
^{------> Module including parameters from the Earth to Mars importation. Will disappear in future dev}

================================================
================================================

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

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

Download in other formats:

Original Format