Context Navigation

inifis_mod.F90 @ 1843

Last change on this file since 1843 was 1843, checked in by jvatant, 7 years ago
Enable tracers management in 1D --JVO
File size: 18.1 KB

Line
1	MODULE inifis_mod
2	IMPLICIT NONE
3
4	CONTAINS
5
6	SUBROUTINE inifis(ngrid,nlayer,nq, &
7	day_ini,pdaysec,nday,ptimestep, &
8	plat,plon,parea, &
9	prad,pg,pr,pcpp)
10
11	use radinc_h, only: ini_radinc_h
12	use datafile_mod
13	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
14	use comgeomfi_h, only: totarea, totarea_planet
15	use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
16	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
17	init_time, daysec, dtphys
18	use comcstfi_mod, only: rad, cpp, g, r, rcp, &
19	mugaz, pi, avocado, kbol
20	use planete_mod, only: nres
21	use planetwide_mod, only: planetwide_sumval
22	use callkeys_mod
23	use mod_phys_lmdz_para, only : is_parallel
24
25	!=======================================================================
26	!
27	! purpose:
28	! -------
29	!
30	! Initialisation for the physical parametrisations of the LMD
31	! Generic Model.
32	!
33	! author: Frederic Hourdin 15 / 10 /93
34	! -------
35	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
36	! Ehouarn Millour (oct. 2008) tracers are now identified
37	! by their names and may not be contiguously
38	! stored in the q(:,:,:,:) array
39	! E.M. (june 2009) use getin routine to load parameters
40	!
41	!
42	! arguments:
43	! ----------
44	!
45	! input:
46	! ------
47	!
48	! ngrid Size of the horizontal grid.
49	! All internal loops are performed on that grid.
50	! nlayer Number of vertical layers.
51	! pdayref Day of reference for the simulation
52	! pday Number of days counted from the North. Spring
53	! equinoxe.
54	!
55	!=======================================================================
56	!
57	!-----------------------------------------------------------------------
58	! declarations:
59	! -------------
60	use ioipsl_getin_p_mod, only: getin_p
61	IMPLICIT NONE
62
63
64
65	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
66	INTEGER,INTENT(IN) :: nday
67	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
68	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
69	integer,intent(in) :: day_ini
70	INTEGER ig,ierr
71
72	EXTERNAL iniorbit,orbite
73	EXTERNAL SSUM
74	REAL SSUM
75
76	! initialize constants in comcstfi_mod
77	rad=prad
78	cpp=pcpp
79	g=pg
80	r=pr
81	rcp=r/cpp
82	pi=2.*asin(1.)
83	avocado = 6.02214179e23 ! added by RW
84	kbol = 1.38064852e-23 ! added by JVO for Titan chem
85
86	! Initialize some "temporal and calendar" related variables
87	CALL init_time(day_ini,pdaysec,nday,ptimestep)
88
89	! read in some parameters from "run.def" for physics,
90	! or shared between dynamics and physics.
91	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
92	! in dynamical steps
93	call getin_p("day_step",day_step) ! number of dynamical steps per day
94	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
95
96	! do we read a startphy.nc file? (default: .true.)
97	call getin_p("startphy_file",startphy_file)
98
99	! --------------------------------------------------------------
100	! Reading the "callphys.def" file controlling some key options
101	! --------------------------------------------------------------
102
103	! check that 'callphys.def' file is around
104	OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
105	CLOSE(99)
106	IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
107
108	!!! IF(ierr.EQ.0) THEN
109	IF(iscallphys) THEN
110	PRINT*
111	PRINT*
112	PRINT*,'--------------------------------------------'
113	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
114	PRINT*,'--------------------------------------------'
115
116	write(,) "Directory where external input files are:"
117	! default 'datadir' is set in "datadir_mod"
118	call getin_p("datadir",datadir) ! default path
119	write(,) " datadir = ",trim(datadir)
120
121	write(,) "Run with or without tracer transport ?"
122	tracer=.false. ! default value
123	call getin_p("tracer",tracer)
124	write(,) " tracer = ",tracer
125
126	write(,) "Run with or without atm mass update ", &
127	" due to tracer evaporation/condensation?"
128	mass_redistrib=.false. ! default value
129	call getin_p("mass_redistrib",mass_redistrib)
130	write(,) " mass_redistrib = ",mass_redistrib
131
132	write(,) "Diurnal cycle ?"
133	write(,) "(if diurnal=false, diurnal averaged solar heating)"
134	diurnal=.true. ! default value
135	call getin_p("diurnal",diurnal)
136	write(,) " diurnal = ",diurnal
137
138	write(,) "Seasonal cycle ?"
139	write(,) "(if season=false, Ls stays constant, to value ", &
140	"set in 'start'"
141	season=.true. ! default value
142	call getin_p("season",season)
143	write(,) " season = ",season
144
145	write(,) "Tidally resonant rotation ?"
146	tlocked=.false. ! default value
147	call getin_p("tlocked",tlocked)
148	write(,) "tlocked = ",tlocked
149
150	write(,) "Saturn ring shadowing ?"
151	rings_shadow = .false.
152	call getin_p("rings_shadow", rings_shadow)
153	write(,) "rings_shadow = ", rings_shadow
154
155	write(,) "Compute latitude-dependent gravity field?"
156	oblate = .false.
157	call getin_p("oblate", oblate)
158	write(,) "oblate = ", oblate
159
160	write(,) "Flattening of the planet (a-b)/a "
161	flatten = 0.0
162	call getin_p("flatten", flatten)
163	write(,) "flatten = ", flatten
164
165	write(,) "Needed if oblate=.true.: J2"
166	J2 = 0.0
167	call getin_p("J2", J2)
168	write(,) "J2 = ", J2
169
170	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
171	MassPlanet = 0.0
172	call getin_p("MassPlanet", MassPlanet)
173	write(,) "MassPlanet = ", MassPlanet
174
175	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
176	Rmean = 0.0
177	call getin_p("Rmean", Rmean)
178	write(,) "Rmean = ", Rmean
179
180	! Test of incompatibility:
181	! if tlocked, then diurnal should be false
182	if (tlocked.and.diurnal) then
183	print*,'If diurnal=true, we should turn off tlocked.'
184	stop
185	endif
186
187	write(,) "Tidal resonance ratio ?"
188	nres=0 ! default value
189	call getin_p("nres",nres)
190	write(,) "nres = ",nres
191
192	write(,) "Write some extra output to the screen ?"
193	lwrite=.false. ! default value
194	call getin_p("lwrite",lwrite)
195	write(,) " lwrite = ",lwrite
196
197	write(,) "Save statistics in file stats.nc ?"
198	callstats=.true. ! default value
199	call getin_p("callstats",callstats)
200	write(,) " callstats = ",callstats
201
202	write(,) "Test energy conservation of model physics ?"
203	enertest=.false. ! default value
204	call getin_p("enertest",enertest)
205	write(,) " enertest = ",enertest
206
207	write(,) "Check to see if cpp values used match gases.def ?"
208	check_cpp_match=.true. ! default value
209	call getin_p("check_cpp_match",check_cpp_match)
210	write(,) " check_cpp_match = ",check_cpp_match
211
212	write(,) "call radiative transfer ?"
213	callrad=.true. ! default value
214	call getin_p("callrad",callrad)
215	write(,) " callrad = ",callrad
216
217	write(,) "call correlated-k radiative transfer ?"
218	corrk=.true. ! default value
219	call getin_p("corrk",corrk)
220	write(,) " corrk = ",corrk
221
222	if (corrk) then
223	! default path is set in datadir
224	write(,) "callcorrk: Correlated-k data base folder:",trim(datadir)
225	call getin_p("corrkdir",corrkdir)
226	write(,) " corrkdir = ",corrkdir
227	endif
228
229	if (corrk .and. ngrid.eq.1) then
230	write(,) "simulate global averaged conditions ?"
231	global1d = .false. ! default value
232	call getin_p("global1d",global1d)
233	write(,) " global1d = ",global1d
234
235	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
236	if (global1d.and.diurnal) then
237	write(,) "if global1d is true, diurnal must be set to false"
238	stop
239	endif
240
241	if (global1d) then
242	write(,) "Solar Zenith angle (deg.) ?"
243	write(,) "(assumed for averaged solar flux S/4)"
244	szangle=60.0 ! default value
245	call getin_p("szangle",szangle)
246	write(,) " szangle = ",szangle
247	endif
248	endif
249
250	write(,) "prohibit calculations outside corrk T grid?"
251	strictboundcorrk=.true. ! default value
252	call getin_p("strictboundcorrk",strictboundcorrk)
253	write(,) "strictboundcorrk = ",strictboundcorrk
254
255	write(,) "call gaseous absorption in the visible bands?", &
256	"(matters only if callrad=T)"
257	callgasvis=.false. ! default value
258	call getin_p("callgasvis",callgasvis)
259	write(,) " callgasvis = ",callgasvis
260
261	write(,) "call continuum opacities in radiative transfer ?", &
262	"(matters only if callrad=T)"
263	continuum=.true. ! default value
264	call getin_p("continuum",continuum)
265	write(,) " continuum = ",continuum
266
267	write(,) "call turbulent vertical diffusion ?"
268	calldifv=.true. ! default value
269	call getin_p("calldifv",calldifv)
270	write(,) " calldifv = ",calldifv
271
272	write(,) "use turbdiff instead of vdifc ?"
273	UseTurbDiff=.true. ! default value
274	call getin_p("UseTurbDiff",UseTurbDiff)
275	write(,) " UseTurbDiff = ",UseTurbDiff
276
277	write(,) "call convective adjustment ?"
278	calladj=.true. ! default value
279	call getin_p("calladj",calladj)
280	write(,) " calladj = ",calladj
281
282	write(,) "Radiative timescale for Newtonian cooling ?"
283	tau_relax=30. ! default value
284	call getin_p("tau_relax",tau_relax)
285	write(,) " tau_relax = ",tau_relax
286	tau_relax=tau_relax243600 ! convert Earth days --> seconds
287
288	write(,)"call thermal conduction in the soil ?"
289	callsoil=.true. ! default value
290	call getin_p("callsoil",callsoil)
291	write(,) " callsoil = ",callsoil
292
293	write(,)"Rad transfer is computed every iradia", &
294	" physical timestep"
295	iradia=1 ! default value
296	call getin_p("iradia",iradia)
297	write(,)" iradia = ",iradia
298
299	write(,)"Rayleigh scattering ?"
300	rayleigh=.false.
301	call getin_p("rayleigh",rayleigh)
302	write(,)" rayleigh = ",rayleigh
303
304	write(,) "Use blackbody for stellar spectrum ?"
305	stelbbody=.false. ! default value
306	call getin_p("stelbbody",stelbbody)
307	write(,) " stelbbody = ",stelbbody
308
309	write(,) "Stellar blackbody temperature ?"
310	stelTbb=5800.0 ! default value
311	call getin_p("stelTbb",stelTbb)
312	write(,) " stelTbb = ",stelTbb
313
314	write(,)"Output mean OLR in 1D?"
315	meanOLR=.false.
316	call getin_p("meanOLR",meanOLR)
317	write(,)" meanOLR = ",meanOLR
318
319	write(,)"Output spectral OLR in 3D?"
320	specOLR=.false.
321	call getin_p("specOLR",specOLR)
322	write(,)" specOLR = ",specOLR
323
324	write(,)"Uniform absorption in radiative transfer?"
325	graybody=.false.
326	call getin_p("graybody",graybody)
327	write(,)" graybody = ",graybody
328
329	! Chemistry
330
331	write(,) "Run with or without chemistry?"
332	callchim=.false. ! default value
333	call getin_p("callchim",callchim)
334	write(,) " callchim = ",callchim
335
336	! sanity check
337	if (callchim.and.(.not.tracer)) then
338	print*,"You are running chemistry without tracer"
339	print*,"Please start again with tracer =.true."
340	stop
341	endif
342
343	write(,)"Chemistry is computed every ichim", &
344	" physical timestep"
345	ichim=1 ! default value
346	call getin_p("ichim",ichim)
347	write(,)" ichim = ",ichim
348
349	! Microphysics
350
351	write(,) "Run with or without microphysics?"
352	callmufi=.false. ! default value
353	call getin_p("callmufi",callmufi)
354	write(,)" callmufi = ",callmufi
355
356	! sanity check
357	if (callmufi.and.(.not.tracer)) then
358	print*,"You are running microphysics without tracer"
359	print*,"Please start again with tracer =.true."
360	stop
361	endif
362
363	write(,) "Compute clouds?"
364	callclouds=.false. ! default value
365	call getin_p("callclouds",callclouds)
366	write(,)" callclouds = ",callclouds
367
368	! sanity check
369	if (callclouds.and.(.not.callmufi)) then
370	print*,"You are trying to make clouds without microphysics !"
371	print*,"Please start again with callmufi =.true."
372	stop
373	endif
374
375	write(,) "Fractal dimension ?"
376	df_mufi=2.0 ! default value
377	call getin_p("df_mufi",df_mufi)
378	write(,)" df_mufi = ",df_mufi
379
380	write(,) "Monomer radius (m) ?"
381	rm_mufi=6.66e-08 ! default value
382	call getin_p("rm_mufi",rm_mufi)
383	write(,)" rm_mufi = ",rm_mufi
384
385	write(,) "Aerosol density (kg.m-3)?"
386	rho_aer_mufi=1.e3 ! default value
387	call getin_p("rho_aer_mufi",rho_aer_mufi)
388	write(,)" rho_aer_mufi = ",rho_aer_mufi
389
390	write(,) "Pressure level of aer. production (Pa) ?"
391	p_prod=1.0 ! default value
392	call getin_p("p_prod",p_prod)
393	write(,)" p_prod = ",p_prod
394
395	write(,) "Aerosol production rate (kg.m-2.s-1) ?"
396	tx_prod=3.5e-13 ! default value
397	call getin_p("tx_prod",tx_prod)
398	write(,)" tx_prod = ",tx_prod
399
400	write(,) "Equivalent radius production (m) ?"
401	rc_prod=2.0e-8 ! default value
402	call getin_p("rc_prod",rc_prod)
403	write(,)" rhc_prod = ",rc_prod
404
405	write(,) "Radius of air (nitrogen) molecule (m) ?"
406	air_rad=1.75e-10 ! default value
407	call getin_p("air_rad",air_rad)
408	write(,)" air_rad = ",air_rad
409
410	write(,) "Path to microphys. config file ?"
411	config_mufi='datagcm/microphysics/config.cfg' ! default value
412	call getin_p("config_mufi",config_mufi)
413	write(,)" config_mufi = ",config_mufi
414
415	! Soil model
416	write(,)"Number of sub-surface layers for soil scheme?"
417	! default value of nsoilmx set in comsoil_h
418	call getin_p("nsoilmx",nsoilmx)
419	write(,)" nsoilmx=",nsoilmx
420
421	write(,)"Thickness of topmost soil layer (m)?"
422	! default value of lay1_soil set in comsoil_h
423	call getin_p("lay1_soil",lay1_soil)
424	write(,)" lay1_soil = ",lay1_soil
425
426	write(,)"Coefficient for soil layer thickness distribution?"
427	! default value of alpha_soil set in comsoil_h
428	call getin_p("alpha_soil",alpha_soil)
429	write(,)" alpha_soil = ",alpha_soil
430
431	write(,)"Remove lower boundary?"
432	nosurf=.false.
433	call getin_p("nosurf",nosurf)
434	write(,)" nosurf = ",nosurf
435
436	! Tests of incompatibility:
437	if (nosurf.and.callsoil) then
438	print*,'nosurf not compatible with soil scheme!'
439	print*,'... got to make a choice!'
440	call abort
441	endif
442
443	write(,)"Add an internal heat flux?", &
444	"... matters only if callsoil=F"
445	intheat=0.
446	call getin_p("intheat",intheat)
447	write(,)" intheat = ",intheat
448
449	write(,)"Use Newtonian cooling for radiative transfer?"
450	newtonian=.false.
451	call getin_p("newtonian",newtonian)
452	write(,)" newtonian = ",newtonian
453
454	! Tests of incompatibility:
455	if (newtonian.and.corrk) then
456	print*,'newtonian not compatible with correlated-k!'
457	call abort
458	endif
459	if (newtonian.and.calladj) then
460	print*,'newtonian not compatible with adjustment!'
461	call abort
462	endif
463	if (newtonian.and.calldifv) then
464	print*,'newtonian not compatible with a boundary layer!'
465	call abort
466	endif
467
468	write(,)"Test physics timescale in 1D?"
469	testradtimes=.false.
470	call getin_p("testradtimes",testradtimes)
471	write(,)" testradtimes = ",testradtimes
472
473	! Test of incompatibility:
474	! if testradtimes used, we must be in 1D
475	if (testradtimes.and.(ngrid.gt.1)) then
476	print*,'testradtimes can only be used in 1D!'
477	call abort
478	endif
479
480	write(,)"Which star?"
481	startype=1 ! default value = Sol
482	call getin_p("startype",startype)
483	write(,)" startype = ",startype
484
485	write(,)"Value of stellar flux at 1 AU?"
486	Fat1AU=1356.0 ! default value = Sol today
487	call getin_p("Fat1AU",Fat1AU)
488	write(,)" Fat1AU = ",Fat1AU
489
490	write(,) "Does user want to force cpp and mugaz?"
491	force_cpp=.false. ! default value
492	call getin_p("force_cpp",force_cpp)
493	write(,) " force_cpp = ",force_cpp
494
495	if (force_cpp) then
496	mugaz = -99999.
497	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
498	call getin_p("mugaz",mugaz)
499	IF (mugaz.eq.-99999.) THEN
500	PRINT *, "mugaz must be set if force_cpp = T"
501	STOP
502	ELSE
503	write(,) "mugaz=",mugaz
504	ENDIF
505	cpp = -99999.
506	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
507	call getin_p("cpp",cpp)
508	IF (cpp.eq.-99999.) THEN
509	PRINT *, "cpp must be set if force_cpp = T"
510	STOP
511	ELSE
512	write(,) "cpp=",cpp
513	ENDIF
514	! else
515	! mugaz=8.314*1000./pr
516	endif ! of if (force_cpp)
517
518
519	call su_gases(nlayer,tracer)
520	call calc_cpp_mugaz
521
522
523	PRINT*,'--------------------------------------------'
524	PRINT*
525	PRINT*
526	ELSE
527	write(,)
528	write(,) 'Cannot read file callphys.def. Is it here ?'
529	stop
530	ENDIF ! of IF(iscallphys)
531
532	PRINT*
533	PRINT*,'inifis: daysec',daysec
534	PRINT*
535	PRINT*,'inifis: The radiative transfer is computed:'
536	PRINT*,' each ',iradia,' physical time-step'
537	PRINT,' or each ',iradiadtphys,' seconds'
538	PRINT*
539	PRINT*,'inifis: Chemistry is computed:'
540	PRINT*,' each ',ichim,' physical time-step'
541	PRINT,' or each ',ichimdtphys,' seconds'
542	PRINT*
543
544	!-----------------------------------------------------------------------
545	! Some more initialization:
546	! ------------------------
547
548	! Initializations for comgeomfi_h
549	totarea=SSUM(ngrid,parea,1)
550	call planetwide_sumval(parea,totarea_planet)
551
552	!! those are defined in comdiurn_h.F90
553	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
554	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
555	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
556	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
557
558	DO ig=1,ngrid
559	sinlat(ig)=sin(plat(ig))
560	coslat(ig)=cos(plat(ig))
561	sinlon(ig)=sin(plon(ig))
562	coslon(ig)=cos(plon(ig))
563	ENDDO
564
565	! initialize variables in radinc_h
566	call ini_radinc_h(nlayer)
567
568	! allocate "comsoil_h" arrays
569	call ini_comsoil_h(ngrid)
570
571	END SUBROUTINE inifis
572
573	END MODULE inifis_mod

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

Download in other formats:

Original Format