Context Navigation

inifis_mod.F90 @ 1862

Last change on this file since 1862 was 1843, checked in by jvatant, 7 years ago
Enable tracers management in 1D --JVO
File size: 18.1 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
[1822]	12	use datafile_mod
[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 ioipsl_getin_p_mod, only: getin_p
	61	IMPLICIT NONE
[1384]	62
[135]	63
	64
[1524]	65	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
[1525]	66	INTEGER,INTENT(IN) :: nday
[1524]	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
[135]	71
[1524]	72	EXTERNAL iniorbit,orbite
	73	EXTERNAL SSUM
	74	REAL SSUM
[135]	75
[1524]	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
[1672]	84	kbol = 1.38064852e-23 ! added by JVO for Titan chem
[135]	85
[1524]	86	! Initialize some "temporal and calendar" related variables
[1525]	87	CALL init_time(day_ini,pdaysec,nday,ptimestep)
[135]	88
[1525]	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
[135]	95
[1670]	96	! do we read a startphy.nc file? (default: .true.)
	97	call getin_p("startphy_file",startphy_file)
	98
[135]	99	! --------------------------------------------------------------
	100	! Reading the "callphys.def" file controlling some key options
	101	! --------------------------------------------------------------
	102
[1524]	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
[135]	107
[1315]	108	!!! IF(ierr.EQ.0) THEN
[1524]	109	IF(iscallphys) THEN
	110	PRINT*
	111	PRINT*
	112	PRINT*,'--------------------------------------------'
	113	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
	114	PRINT*,'--------------------------------------------'
[135]	115
[1524]	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)
[374]	120
[1524]	121	write(,) "Run with or without tracer transport ?"
	122	tracer=.false. ! default value
	123	call getin_p("tracer",tracer)
	124	write(,) " tracer = ",tracer
[135]	125
[1524]	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
[728]	131
[1524]	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
[135]	137
[1524]	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
[135]	144
[1524]	145	write(,) "Tidally resonant rotation ?"
	146	tlocked=.false. ! default value
	147	call getin_p("tlocked",tlocked)
	148	write(,) "tlocked = ",tlocked
[135]	149
[1524]	150	write(,) "Saturn ring shadowing ?"
	151	rings_shadow = .false.
	152	call getin_p("rings_shadow", rings_shadow)
	153	write(,) "rings_shadow = ", rings_shadow
[1133]	154
[1524]	155	write(,) "Compute latitude-dependent gravity field?"
	156	oblate = .false.
	157	call getin_p("oblate", oblate)
	158	write(,) "oblate = ", oblate
[1194]	159
[1524]	160	write(,) "Flattening of the planet (a-b)/a "
	161	flatten = 0.0
	162	call getin_p("flatten", flatten)
[1672]	163	write(,) "flatten = ", flatten
[1194]	164
[1524]	165	write(,) "Needed if oblate=.true.: J2"
	166	J2 = 0.0
	167	call getin_p("J2", J2)
	168	write(,) "J2 = ", J2
[1194]	169
[1524]	170	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
	171	MassPlanet = 0.0
	172	call getin_p("MassPlanet", MassPlanet)
	173	write(,) "MassPlanet = ", MassPlanet
[1194]	174
[1524]	175	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
	176	Rmean = 0.0
	177	call getin_p("Rmean", Rmean)
	178	write(,) "Rmean = ", Rmean
[1194]	179
[135]	180	! Test of incompatibility:
	181	! if tlocked, then diurnal should be false
[1524]	182	if (tlocked.and.diurnal) then
	183	print*,'If diurnal=true, we should turn off tlocked.'
	184	stop
	185	endif
[135]	186
[1524]	187	write(,) "Tidal resonance ratio ?"
	188	nres=0 ! default value
	189	call getin_p("nres",nres)
	190	write(,) "nres = ",nres
[135]	191
[1524]	192	write(,) "Write some extra output to the screen ?"
	193	lwrite=.false. ! default value
	194	call getin_p("lwrite",lwrite)
	195	write(,) " lwrite = ",lwrite
[135]	196
[1524]	197	write(,) "Save statistics in file stats.nc ?"
	198	callstats=.true. ! default value
	199	call getin_p("callstats",callstats)
	200	write(,) " callstats = ",callstats
[135]	201
[1524]	202	write(,) "Test energy conservation of model physics ?"
	203	enertest=.false. ! default value
	204	call getin_p("enertest",enertest)
	205	write(,) " enertest = ",enertest
[135]	206
[1524]	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
[538]	211
[1524]	212	write(,) "call radiative transfer ?"
	213	callrad=.true. ! default value
	214	call getin_p("callrad",callrad)
	215	write(,) " callrad = ",callrad
[135]	216
[1524]	217	write(,) "call correlated-k radiative transfer ?"
	218	corrk=.true. ! default value
	219	call getin_p("corrk",corrk)
	220	write(,) " corrk = ",corrk
[1822]	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
[135]	240
[1822]	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
[1524]	250	write(,) "prohibit calculations outside corrk T grid?"
	251	strictboundcorrk=.true. ! default value
	252	call getin_p("strictboundcorrk",strictboundcorrk)
	253	write(,) "strictboundcorrk = ",strictboundcorrk
[1145]	254
[1524]	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
[538]	260
[1524]	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
[538]	266
[1524]	267	write(,) "call turbulent vertical diffusion ?"
	268	calldifv=.true. ! default value
	269	call getin_p("calldifv",calldifv)
	270	write(,) " calldifv = ",calldifv
[135]	271
[1524]	272	write(,) "use turbdiff instead of vdifc ?"
	273	UseTurbDiff=.true. ! default value
	274	call getin_p("UseTurbDiff",UseTurbDiff)
	275	write(,) " UseTurbDiff = ",UseTurbDiff
[596]	276
[1524]	277	write(,) "call convective adjustment ?"
	278	calladj=.true. ! default value
	279	call getin_p("calladj",calladj)
	280	write(,) " calladj = ",calladj
[135]	281
[1524]	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
[253]	287
[1524]	288	write(,)"call thermal conduction in the soil ?"
	289	callsoil=.true. ! default value
	290	call getin_p("callsoil",callsoil)
	291	write(,) " callsoil = ",callsoil
[135]	292
[1524]	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
[135]	298
[1524]	299	write(,)"Rayleigh scattering ?"
	300	rayleigh=.false.
	301	call getin_p("rayleigh",rayleigh)
	302	write(,)" rayleigh = ",rayleigh
[135]	303
[1524]	304	write(,) "Use blackbody for stellar spectrum ?"
	305	stelbbody=.false. ! default value
	306	call getin_p("stelbbody",stelbbody)
	307	write(,) " stelbbody = ",stelbbody
[135]	308
[1524]	309	write(,) "Stellar blackbody temperature ?"
	310	stelTbb=5800.0 ! default value
	311	call getin_p("stelTbb",stelTbb)
	312	write(,) " stelTbb = ",stelTbb
[253]	313
[1524]	314	write(,)"Output mean OLR in 1D?"
	315	meanOLR=.false.
	316	call getin_p("meanOLR",meanOLR)
	317	write(,)" meanOLR = ",meanOLR
[135]	318
[1524]	319	write(,)"Output spectral OLR in 3D?"
	320	specOLR=.false.
	321	call getin_p("specOLR",specOLR)
	322	write(,)" specOLR = ",specOLR
[135]	323
[1524]	324	write(,)"Uniform absorption in radiative transfer?"
	325	graybody=.false.
	326	call getin_p("graybody",graybody)
	327	write(,)" graybody = ",graybody
[253]	328
[1672]	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
[1795]	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)
[1822]	378	write(,)" df_mufi = ",df_mufi
[1795]	379
	380	write(,) "Monomer radius (m) ?"
	381	rm_mufi=6.66e-08 ! default value
	382	call getin_p("rm_mufi",rm_mufi)
[1822]	383	write(,)" rm_mufi = ",rm_mufi
[1795]	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)
[1822]	388	write(,)" rho_aer_mufi = ",rho_aer_mufi
[1795]	389
	390	write(,) "Pressure level of aer. production (Pa) ?"
	391	p_prod=1.0 ! default value
	392	call getin_p("p_prod",p_prod)
[1822]	393	write(,)" p_prod = ",p_prod
[1795]	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)
[1822]	398	write(,)" tx_prod = ",tx_prod
[1795]	399
	400	write(,) "Equivalent radius production (m) ?"
	401	rc_prod=2.0e-8 ! default value
	402	call getin_p("rc_prod",rc_prod)
[1822]	403	write(,)" rhc_prod = ",rc_prod
[1795]	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)
[1822]	408	write(,)" air_rad = ",air_rad
[1795]	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)
[1822]	413	write(,)" config_mufi = ",config_mufi
[1795]	414
[1538]	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
[1524]	431	write(,)"Remove lower boundary?"
	432	nosurf=.false.
	433	call getin_p("nosurf",nosurf)
	434	write(,)" nosurf = ",nosurf
[253]	435
	436	! Tests of incompatibility:
[1524]	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
[253]	442
[1524]	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
[952]	448
[1524]	449	write(,)"Use Newtonian cooling for radiative transfer?"
	450	newtonian=.false.
	451	call getin_p("newtonian",newtonian)
	452	write(,)" newtonian = ",newtonian
[253]	453
	454	! Tests of incompatibility:
[1524]	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
[253]	467
[1524]	468	write(,)"Test physics timescale in 1D?"
	469	testradtimes=.false.
	470	call getin_p("testradtimes",testradtimes)
	471	write(,)" testradtimes = ",testradtimes
[253]	472
	473	! Test of incompatibility:
	474	! if testradtimes used, we must be in 1D
[1524]	475	if (testradtimes.and.(ngrid.gt.1)) then
	476	print*,'testradtimes can only be used in 1D!'
	477	call abort
	478	endif
[253]	479
[1524]	480	write(,)"Which star?"
	481	startype=1 ! default value = Sol
	482	call getin_p("startype",startype)
	483	write(,)" startype = ",startype
[135]	484
[1524]	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
[135]	489
[1524]	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
[589]	494
[1524]	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
[961]	514	! else
	515	! mugaz=8.314*1000./pr
[1524]	516	endif ! of if (force_cpp)
[1648]	517
	518
	519	call su_gases(nlayer,tracer)
[1524]	520	call calc_cpp_mugaz
[1648]	521
	522
[1524]	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)
[135]	531
[1524]	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*
[1672]	539	PRINT*,'inifis: Chemistry is computed:'
	540	PRINT*,' each ',ichim,' physical time-step'
	541	PRINT,' or each ',ichimdtphys,' seconds'
	542	PRINT*
[135]	543
	544	!-----------------------------------------------------------------------
	545	! Some more initialization:
	546	! ------------------------
	547
[1542]	548	! Initializations for comgeomfi_h
	549	totarea=SSUM(ngrid,parea,1)
	550	call planetwide_sumval(parea,totarea_planet)
[787]	551
[1524]	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))
[787]	557
[1524]	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
[1216]	564
[1722]	565	! initialize variables in radinc_h
[1529]	566	call ini_radinc_h(nlayer)
	567
[1524]	568	! allocate "comsoil_h" arrays
	569	call ini_comsoil_h(ngrid)
	570
	571	END SUBROUTINE inifis
[135]	572
[1524]	573	END MODULE inifis_mod

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

Download in other formats:

Original Format