Context Navigation

inifis_mod.F90 @ 1819

Last change on this file since 1819 was 1795, checked in by jvatant, 8 years ago
Making Titan's hazy again, part II + Added calmufi and inimufi routines that interface YAMMS model + Major changes of the tracer gestion in tracer_h (new query by name) + Update the deftank JVO
File size: 16.9 KB

Rev	Line
[1524]	1	MODULE inifis_mod
	2	IMPLICIT NONE
[253]	3
[1524]	4	CONTAINS
[253]	5
[1524]	6	SUBROUTINE inifis(ngrid,nlayer,nq, &
[1525]	7	day_ini,pdaysec,nday,ptimestep, &
[1524]	8	plat,plon,parea, &
	9	prad,pg,pr,pcpp)
	10
[1788]	11	use radinc_h, only: ini_radinc_h
[1795]	12	use datafile_mod, only: datadir,config_mufi
[1524]	13	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
[1542]	14	use comgeomfi_h, only: totarea, totarea_planet
[1538]	15	use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
[1525]	16	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
	17	init_time, daysec, dtphys
[1524]	18	use comcstfi_mod, only: rad, cpp, g, r, rcp, &
[1672]	19	mugaz, pi, avocado, kbol
[1524]	20	use planete_mod, only: nres
	21	use planetwide_mod, only: planetwide_sumval
	22	use callkeys_mod
[1529]	23	use mod_phys_lmdz_para, only : is_parallel
[1524]	24
[135]	25	!=======================================================================
	26	!
	27	! purpose:
	28	! -------
	29	!
	30	! Initialisation for the physical parametrisations of the LMD
[305]	31	! Generic Model.
[135]	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	! -------------
[1524]	60	use datafile_mod, only: datadir
	61	use ioipsl_getin_p_mod, only: getin_p
	62	IMPLICIT NONE
[1384]	63
[135]	64
	65
[1524]	66	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
[1525]	67	INTEGER,INTENT(IN) :: nday
[1524]	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
[135]	72
[1524]	73	EXTERNAL iniorbit,orbite
	74	EXTERNAL SSUM
	75	REAL SSUM
[135]	76
[1524]	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
[1672]	85	kbol = 1.38064852e-23 ! added by JVO for Titan chem
[135]	86
[1524]	87	! Initialize some "temporal and calendar" related variables
[1525]	88	CALL init_time(day_ini,pdaysec,nday,ptimestep)
[135]	89
[1525]	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
[135]	96
[1670]	97	! do we read a startphy.nc file? (default: .true.)
	98	call getin_p("startphy_file",startphy_file)
	99
[135]	100	! --------------------------------------------------------------
	101	! Reading the "callphys.def" file controlling some key options
	102	! --------------------------------------------------------------
	103
[1524]	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
[135]	108
[1315]	109	!!! IF(ierr.EQ.0) THEN
[1524]	110	IF(iscallphys) THEN
	111	PRINT*
	112	PRINT*
	113	PRINT*,'--------------------------------------------'
	114	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
	115	PRINT*,'--------------------------------------------'
[135]	116
[1524]	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)
[374]	121
[1524]	122	write(,) "Run with or without tracer transport ?"
	123	tracer=.false. ! default value
	124	call getin_p("tracer",tracer)
	125	write(,) " tracer = ",tracer
[135]	126
[1524]	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
[728]	132
[1524]	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
[135]	138
[1524]	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
[135]	145
[1524]	146	write(,) "Tidally resonant rotation ?"
	147	tlocked=.false. ! default value
	148	call getin_p("tlocked",tlocked)
	149	write(,) "tlocked = ",tlocked
[135]	150
[1524]	151	write(,) "Saturn ring shadowing ?"
	152	rings_shadow = .false.
	153	call getin_p("rings_shadow", rings_shadow)
	154	write(,) "rings_shadow = ", rings_shadow
[1133]	155
[1524]	156	write(,) "Compute latitude-dependent gravity field?"
	157	oblate = .false.
	158	call getin_p("oblate", oblate)
	159	write(,) "oblate = ", oblate
[1194]	160
[1524]	161	write(,) "Flattening of the planet (a-b)/a "
	162	flatten = 0.0
	163	call getin_p("flatten", flatten)
[1672]	164	write(,) "flatten = ", flatten
[1194]	165
[1524]	166	write(,) "Needed if oblate=.true.: J2"
	167	J2 = 0.0
	168	call getin_p("J2", J2)
	169	write(,) "J2 = ", J2
[1194]	170
[1524]	171	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
	172	MassPlanet = 0.0
	173	call getin_p("MassPlanet", MassPlanet)
	174	write(,) "MassPlanet = ", MassPlanet
[1194]	175
[1524]	176	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
	177	Rmean = 0.0
	178	call getin_p("Rmean", Rmean)
	179	write(,) "Rmean = ", Rmean
[1194]	180
[135]	181	! Test of incompatibility:
	182	! if tlocked, then diurnal should be false
[1524]	183	if (tlocked.and.diurnal) then
	184	print*,'If diurnal=true, we should turn off tlocked.'
	185	stop
	186	endif
[135]	187
[1524]	188	write(,) "Tidal resonance ratio ?"
	189	nres=0 ! default value
	190	call getin_p("nres",nres)
	191	write(,) "nres = ",nres
[135]	192
[1524]	193	write(,) "Write some extra output to the screen ?"
	194	lwrite=.false. ! default value
	195	call getin_p("lwrite",lwrite)
	196	write(,) " lwrite = ",lwrite
[135]	197
[1524]	198	write(,) "Save statistics in file stats.nc ?"
	199	callstats=.true. ! default value
	200	call getin_p("callstats",callstats)
	201	write(,) " callstats = ",callstats
[135]	202
[1524]	203	write(,) "Test energy conservation of model physics ?"
	204	enertest=.false. ! default value
	205	call getin_p("enertest",enertest)
	206	write(,) " enertest = ",enertest
[135]	207
[1524]	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
[538]	212
[1524]	213	write(,) "call radiative transfer ?"
	214	callrad=.true. ! default value
	215	call getin_p("callrad",callrad)
	216	write(,) " callrad = ",callrad
[135]	217
[1524]	218	write(,) "call correlated-k radiative transfer ?"
	219	corrk=.true. ! default value
	220	call getin_p("corrk",corrk)
	221	write(,) " corrk = ",corrk
[135]	222
[1524]	223	write(,) "prohibit calculations outside corrk T grid?"
	224	strictboundcorrk=.true. ! default value
	225	call getin_p("strictboundcorrk",strictboundcorrk)
	226	write(,) "strictboundcorrk = ",strictboundcorrk
[1145]	227
[1524]	228	write(,) "call gaseous absorption in the visible bands?", &
	229	"(matters only if callrad=T)"
	230	callgasvis=.false. ! default value
	231	call getin_p("callgasvis",callgasvis)
	232	write(,) " callgasvis = ",callgasvis
[538]	233
[1524]	234	write(,) "call continuum opacities in radiative transfer ?", &
	235	"(matters only if callrad=T)"
	236	continuum=.true. ! default value
	237	call getin_p("continuum",continuum)
	238	write(,) " continuum = ",continuum
[538]	239
[1524]	240	write(,) "call turbulent vertical diffusion ?"
	241	calldifv=.true. ! default value
	242	call getin_p("calldifv",calldifv)
	243	write(,) " calldifv = ",calldifv
[135]	244
[1524]	245	write(,) "use turbdiff instead of vdifc ?"
	246	UseTurbDiff=.true. ! default value
	247	call getin_p("UseTurbDiff",UseTurbDiff)
	248	write(,) " UseTurbDiff = ",UseTurbDiff
[596]	249
[1524]	250	write(,) "call convective adjustment ?"
	251	calladj=.true. ! default value
	252	call getin_p("calladj",calladj)
	253	write(,) " calladj = ",calladj
[135]	254
[1524]	255	write(,) "Radiative timescale for Newtonian cooling ?"
	256	tau_relax=30. ! default value
	257	call getin_p("tau_relax",tau_relax)
	258	write(,) " tau_relax = ",tau_relax
	259	tau_relax=tau_relax243600 ! convert Earth days --> seconds
[253]	260
[1524]	261	write(,)"call thermal conduction in the soil ?"
	262	callsoil=.true. ! default value
	263	call getin_p("callsoil",callsoil)
	264	write(,) " callsoil = ",callsoil
[135]	265
[1524]	266	write(,)"Rad transfer is computed every iradia", &
	267	" physical timestep"
	268	iradia=1 ! default value
	269	call getin_p("iradia",iradia)
	270	write(,)" iradia = ",iradia
[135]	271
[1524]	272	write(,)"Rayleigh scattering ?"
	273	rayleigh=.false.
	274	call getin_p("rayleigh",rayleigh)
	275	write(,)" rayleigh = ",rayleigh
[135]	276
[1524]	277	write(,) "Use blackbody for stellar spectrum ?"
	278	stelbbody=.false. ! default value
	279	call getin_p("stelbbody",stelbbody)
	280	write(,) " stelbbody = ",stelbbody
[135]	281
[1524]	282	write(,) "Stellar blackbody temperature ?"
	283	stelTbb=5800.0 ! default value
	284	call getin_p("stelTbb",stelTbb)
	285	write(,) " stelTbb = ",stelTbb
[253]	286
[1524]	287	write(,)"Output mean OLR in 1D?"
	288	meanOLR=.false.
	289	call getin_p("meanOLR",meanOLR)
	290	write(,)" meanOLR = ",meanOLR
[135]	291
[1524]	292	write(,)"Output spectral OLR in 3D?"
	293	specOLR=.false.
	294	call getin_p("specOLR",specOLR)
	295	write(,)" specOLR = ",specOLR
[135]	296
[1524]	297	write(,)"Uniform absorption in radiative transfer?"
	298	graybody=.false.
	299	call getin_p("graybody",graybody)
	300	write(,)" graybody = ",graybody
[253]	301
[1672]	302	! Chemistry
	303
	304	write(,) "Run with or without chemistry?"
	305	callchim=.false. ! default value
	306	call getin_p("callchim",callchim)
	307	write(,) " callchim = ",callchim
	308
	309	! sanity check
	310	if (callchim.and.(.not.tracer)) then
	311	print*,"You are running chemistry without tracer"
	312	print*,"Please start again with tracer =.true."
	313	stop
	314	endif
	315
	316	write(,)"Chemistry is computed every ichim", &
	317	" physical timestep"
	318	ichim=1 ! default value
	319	call getin_p("ichim",ichim)
	320	write(,)" ichim = ",ichim
	321
[1795]	322	! Microphysics
	323
	324	write(,) "Run with or without microphysics?"
	325	callmufi=.false. ! default value
	326	call getin_p("callmufi",callmufi)
	327	write(,)" callmufi = ",callmufi
	328
	329	! sanity check
	330	if (callmufi.and.(.not.tracer)) then
	331	print*,"You are running microphysics without tracer"
	332	print*,"Please start again with tracer =.true."
	333	stop
	334	endif
	335
	336	write(,) "Compute clouds?"
	337	callclouds=.false. ! default value
	338	call getin_p("callclouds",callclouds)
	339	write(,)" callclouds = ",callclouds
	340
	341	! sanity check
	342	if (callclouds.and.(.not.callmufi)) then
	343	print*,"You are trying to make clouds without microphysics !"
	344	print*,"Please start again with callmufi =.true."
	345	stop
	346	endif
	347
	348	write(,) "Fractal dimension ?"
	349	df_mufi=2.0 ! default value
	350	call getin_p("df_mufi",df_mufi)
	351
	352	write(,) "Monomer radius (m) ?"
	353	rm_mufi=6.66e-08 ! default value
	354	call getin_p("rm_mufi",rm_mufi)
	355
	356	write(,) "Aerosol density (kg.m-3)?"
	357	rho_aer_mufi=1.e3 ! default value
	358	call getin_p("rho_aer_mufi",rho_aer_mufi)
	359
	360	write(,) "Pressure level of aer. production (Pa) ?"
	361	p_prod=1.0 ! default value
	362	call getin_p("p_prod",p_prod)
	363
	364	write(,) "Aerosol production rate (kg.m-2.s-1) ?"
	365	tx_prod=3.5e-13 ! default value
	366	call getin_p("tx_prod",tx_prod)
	367
	368	write(,) "Equivalent radius production (m) ?"
	369	rc_prod=2.0e-8 ! default value
	370	call getin_p("rc_prod",rc_prod)
	371
	372	write(,) "Radius of air (nitrogen) molecule (m) ?"
	373	air_rad=1.75e-10 ! default value
	374	call getin_p("air_rad",air_rad)
	375
	376	write(,) "Path to microphys. config file ?"
	377	config_mufi='datagcm/microphysics/config.cfg' ! default value
	378	call getin_p("config_mufi",config_mufi)
	379
[1538]	380	! Soil model
	381	write(,)"Number of sub-surface layers for soil scheme?"
	382	! default value of nsoilmx set in comsoil_h
	383	call getin_p("nsoilmx",nsoilmx)
	384	write(,)" nsoilmx=",nsoilmx
	385
	386	write(,)"Thickness of topmost soil layer (m)?"
	387	! default value of lay1_soil set in comsoil_h
	388	call getin_p("lay1_soil",lay1_soil)
	389	write(,)" lay1_soil = ",lay1_soil
	390
	391	write(,)"Coefficient for soil layer thickness distribution?"
	392	! default value of alpha_soil set in comsoil_h
	393	call getin_p("alpha_soil",alpha_soil)
	394	write(,)" alpha_soil = ",alpha_soil
	395
[1524]	396	write(,)"Remove lower boundary?"
	397	nosurf=.false.
	398	call getin_p("nosurf",nosurf)
	399	write(,)" nosurf = ",nosurf
[253]	400
	401	! Tests of incompatibility:
[1524]	402	if (nosurf.and.callsoil) then
	403	print*,'nosurf not compatible with soil scheme!'
	404	print*,'... got to make a choice!'
	405	call abort
	406	endif
[253]	407
[1524]	408	write(,)"Add an internal heat flux?", &
	409	"... matters only if callsoil=F"
	410	intheat=0.
	411	call getin_p("intheat",intheat)
	412	write(,)" intheat = ",intheat
[952]	413
[1524]	414	write(,)"Use Newtonian cooling for radiative transfer?"
	415	newtonian=.false.
	416	call getin_p("newtonian",newtonian)
	417	write(,)" newtonian = ",newtonian
[253]	418
	419	! Tests of incompatibility:
[1524]	420	if (newtonian.and.corrk) then
	421	print*,'newtonian not compatible with correlated-k!'
	422	call abort
	423	endif
	424	if (newtonian.and.calladj) then
	425	print*,'newtonian not compatible with adjustment!'
	426	call abort
	427	endif
	428	if (newtonian.and.calldifv) then
	429	print*,'newtonian not compatible with a boundary layer!'
	430	call abort
	431	endif
[253]	432
[1524]	433	write(,)"Test physics timescale in 1D?"
	434	testradtimes=.false.
	435	call getin_p("testradtimes",testradtimes)
	436	write(,)" testradtimes = ",testradtimes
[253]	437
	438	! Test of incompatibility:
	439	! if testradtimes used, we must be in 1D
[1524]	440	if (testradtimes.and.(ngrid.gt.1)) then
	441	print*,'testradtimes can only be used in 1D!'
	442	call abort
	443	endif
[253]	444
[1524]	445	write(,)"Which star?"
	446	startype=1 ! default value = Sol
	447	call getin_p("startype",startype)
	448	write(,)" startype = ",startype
[135]	449
[1524]	450	write(,)"Value of stellar flux at 1 AU?"
	451	Fat1AU=1356.0 ! default value = Sol today
	452	call getin_p("Fat1AU",Fat1AU)
	453	write(,)" Fat1AU = ",Fat1AU
[135]	454
[1524]	455	write(,) "Does user want to force cpp and mugaz?"
	456	force_cpp=.false. ! default value
	457	call getin_p("force_cpp",force_cpp)
	458	write(,) " force_cpp = ",force_cpp
[589]	459
[1524]	460	if (force_cpp) then
	461	mugaz = -99999.
	462	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
	463	call getin_p("mugaz",mugaz)
	464	IF (mugaz.eq.-99999.) THEN
	465	PRINT *, "mugaz must be set if force_cpp = T"
	466	STOP
	467	ELSE
	468	write(,) "mugaz=",mugaz
	469	ENDIF
	470	cpp = -99999.
	471	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
	472	call getin_p("cpp",cpp)
	473	IF (cpp.eq.-99999.) THEN
	474	PRINT *, "cpp must be set if force_cpp = T"
	475	STOP
	476	ELSE
	477	write(,) "cpp=",cpp
	478	ENDIF
[961]	479	! else
	480	! mugaz=8.314*1000./pr
[1524]	481	endif ! of if (force_cpp)
[1648]	482
	483
	484	call su_gases(nlayer,tracer)
[1524]	485	call calc_cpp_mugaz
[1648]	486
	487
[1524]	488	PRINT*,'--------------------------------------------'
	489	PRINT*
	490	PRINT*
	491	ELSE
	492	write(,)
	493	write(,) 'Cannot read file callphys.def. Is it here ?'
	494	stop
	495	ENDIF ! of IF(iscallphys)
[135]	496
[1524]	497	PRINT*
	498	PRINT*,'inifis: daysec',daysec
	499	PRINT*
	500	PRINT*,'inifis: The radiative transfer is computed:'
	501	PRINT*,' each ',iradia,' physical time-step'
	502	PRINT,' or each ',iradiadtphys,' seconds'
	503	PRINT*
[1672]	504	PRINT*,'inifis: Chemistry is computed:'
	505	PRINT*,' each ',ichim,' physical time-step'
	506	PRINT,' or each ',ichimdtphys,' seconds'
	507	PRINT*
[135]	508
	509	!-----------------------------------------------------------------------
	510	! Some more initialization:
	511	! ------------------------
	512
[1542]	513	! Initializations for comgeomfi_h
	514	totarea=SSUM(ngrid,parea,1)
	515	call planetwide_sumval(parea,totarea_planet)
[787]	516
[1524]	517	!! those are defined in comdiurn_h.F90
	518	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
	519	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
	520	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
	521	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
[787]	522
[1524]	523	DO ig=1,ngrid
	524	sinlat(ig)=sin(plat(ig))
	525	coslat(ig)=cos(plat(ig))
	526	sinlon(ig)=sin(plon(ig))
	527	coslon(ig)=cos(plon(ig))
	528	ENDDO
[1216]	529
[1722]	530	! initialize variables in radinc_h
[1529]	531	call ini_radinc_h(nlayer)
	532
[1524]	533	! allocate "comsoil_h" arrays
	534	call ini_comsoil_h(ngrid)
	535
	536	END SUBROUTINE inifis
[135]	537
[1524]	538	END MODULE inifis_mod

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

Download in other formats:

Original Format