Context Navigation

inifis_mod.F90 @ 1822

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

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

Download in other formats:

Original Format