Context Navigation

tabfi.F @ 3026

Last change on this file since 3026 was 2311, checked in by emillour, 5 years ago
Mars GCM: Code tidying: use getin_p() instead of getin() and use "call abort_physic" instead of "stop" or "call abort". EM
File size: 23.1 KB

Rev	Line
[38]	1	c=======================================================================
	2	SUBROUTINE tabfi(nid,Lmodif,tab0,day_ini,lmax,p_rad,
	3	. p_omeg,p_g,p_mugaz,p_daysec,time)
	4	c=======================================================================
	5	c
	6	c C. Hourdin 15/11/96
	7	c
	8	c Object: Lecture du tab_cntrl physique dans un fichier
	9	c ------ et initialisation des constantes physiques
	10	c
	11	c Arguments:
	12	c ----------
	13	c
	14	c Inputs:
	15	c ------
	16	c
	17	c - nid: unitne logique du fichier ou on va lire le tab_cntrl
	18	c (ouvert dans le programme appellant)
	19	c
	20	c si nid=0:
	21	c pas de lecture du tab_cntrl mais
	22	c Valeurs par default des constantes physiques
	23	c
	24	c - tab0: Offset de tab_cntrl a partir duquel sont ranges
	25	c les parametres physiques (50 pour start_archive)
	26	c
	27	c - Lmodif: si on souhaite modifier les constantes Lmodif = 1 = TRUE
	28	c
	29	c
	30	c Outputs:
	31	c --------
	32	c
	33	c - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite
	34	c comparer avec le day_ini dynamique)
	35	c
[1266]	36	c - lmax: tab_cntrl(tab0+2) (pour test avec nlayer)
[38]	37	c
	38	c - p_rad
	39	c - p_omeg !
	40	c - p_g ! Constantes physiques ayant des
	41	c - p_mugaz ! homonymes dynamiques
	42	c - p_daysec !
	43	c
	44	c=======================================================================
[2281]	45	! to use 'getin_p'
	46	use ioipsl_getin_p_mod, only: getin_p
[1047]	47
	48	use comsoil_h, only: volcapa ! soil volumetric heat capacity
	49	use surfdat_h, only: z0_default, emissiv, emisice, albedice,
	50	& iceradius, dtemisice, iceradius
[1246]	51	use dimradmars_mod, only: tauvis
[1130]	52	use iostart, only: get_var
	53	use mod_phys_lmdz_para, only: is_parallel
[1524]	54	use comcstfi_h, only: g, mugaz, omeg, rad, rcp
	55	use time_phylmdz_mod, only: daysec, dtphys
	56	use planete_h, only: aphelie, emin_turb, lmixmin, obliquit,
	57	& peri_day, periheli, year_day
[38]	58	implicit none
	59
[1543]	60	include "netcdf.inc"
[38]	61
	62	c-----------------------------------------------------------------------
	63	c Declarations
	64	c-----------------------------------------------------------------------
	65
	66	c Arguments
	67	c ---------
[1130]	68	INTEGER,INTENT(IN) :: nid,tab0
	69	INTEGER*4,INTENT(OUT) :: day_ini
	70	INTEGER,INTENT(IN) :: Lmodif
	71	INTEGER,INTENT(OUT) :: lmax
	72	REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_mugaz,p_daysec,time
[38]	73
	74	c Variables
	75	c ---------
[1130]	76	INTEGER :: nvarid
	77	REAL :: peri_ls
[38]	78	INTEGER length
	79	parameter (length = 100)
	80	REAL tab_cntrl(length) ! array in which are stored the run's parameters
	81	INTEGER ierr
	82	INTEGER size
	83	CHARACTER modif*20
[1130]	84	LOGICAL :: found
[2281]	85	CHARACTER(len=5) :: modname="tabfi"
[38]	86
[1130]	87	write(,)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
	88
[38]	89	c-----------------------------------------------------------------------
	90	c Initialization of various physical constants to defaut values (nid = 0 case)
	91	c-----------------------------------------------------------------------
	92	IF (nid.eq.0) then
[2281]	93
	94	! to avoid further issues with writing
	95	tab_cntrl(:)=0
	96	lmax=0
	97	day_ini=0
	98	time = 0
[38]	99
	100	c Reference pressure
	101	c-------------------------------------
	102	c pressrf = 670. ! Pression de reference (Pa) ~650
	103
	104	c Infos about Mars for the dynamics and physics
	105	c----------------------------------------------------------
	106	rad=3397200. ! radius of Mars (m) ~3397200 m
	107	daysec=88775. ! length of a sol (s) ~88775 s
	108	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
	109	g=3.72 ! gravity (m.s-2) ~3.72
	110	mugaz=43.49 ! Molar mass of the atmosphere (g.mol-1) ~43.49
	111	rcp=.256793 ! = r/cp ~0.256793
	112
	113	c Informations about Mars, only for physics
	114	c-----------------------------------------------------
	115	year_day = 669. !Modif FH: length of year (sols) ~668.6
	116	periheli = 206.66 ! min. Sun-Mars distance (Mkm) ~206.66
	117	aphelie = 249.22 ! max. Sun-Mars distance (Mkm) ~249.22
	118	peri_day = 485. ! date of perihelion (sols since N. spring)
	119	obliquit = 25.19 ! Obliquity of the planet (deg) ~25.19
	120
	121	c Boundary layer and turbulence
	122	c----------------------------
[224]	123	z0_default = 1.e-2 ! surface roughness (m) ~0.01
[38]	124	emin_turb = 1.e-6 ! minimal energy ~1.e-8
	125	lmixmin = 30 ! mixing length ~100
	126
	127	c Optical properties of polar caps and ground emissivity
	128	c-----------------------------------------------------
	129	emissiv=.95 ! Emissivity of martian soil ~.95
	130	emisice(1)=0.95 ! Emissivity of northern cap
	131	emisice(2)=0.95 ! Emissivity of southern cap
	132	albedice(1)=0.65 ! Albedo of northern cap
	133	albedice(2)=0.65 ! Albedo of southern cap
	134	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	135	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	136	dtemisice(1) = 0.4 ! time scale for snow metamorphism (north)
	137	dtemisice(2) = 0.4 ! time scale for snow metamorphism (south)
	138
	139	c dust aerosol properties
	140	c---------------------------------
	141	tauvis= 0.2 ! mean visible optical depth
	142
	143	c Ancien code radiatif (non utilise avec le code d'apres 03/96)
	144	c---------------------------------------------------------------
	145	c tauir= 0. ! .2 ratio (mean IR opt.depth)/Visible
	146	c scatalb=0. ! .86 scaterring albedo visible (~.86)
	147	c asfact=0. ! .79 assymetrie factor visible (~.79)
	148	c day0 = 0 ! = 0 en general !!!
	149
	150	c soil properties
	151	volcapa = 1.e6 ! soil volumetric heat capacity (in comsoil.h)
	152	ELSE
	153	c-----------------------------------------------------------------------
	154	c Initialization of physical constants by reading array tab_cntrl(:)
	155	c which contains these parameters (nid != 0 case)
	156	c-----------------------------------------------------------------------
	157	c Read 'controle' array
	158	c
[1130]	159	call get_var("controle",tab_cntrl,found)
	160	if (.not.found) then
[2281]	161	call abort_physic(modname,
	162	& "tabfi: Failed reading <controle> array",1)
[1130]	163	else
	164	write(,)'tabfi: tab_cntrl',tab_cntrl
	165	endif
[38]	166	c
	167	c Initialization of some physical constants
	168	c informations on physics grid
	169	lmax = nint(tab_cntrl(tab0+2))
	170	day_ini = tab_cntrl(tab0+3)
	171	time = tab_cntrl(tab0+4)
	172	write (,) 'IN tabfi day_ini=',day_ini
	173	c Informations about planet Mars for dynamics and physics
	174	rad = tab_cntrl(tab0+5)
	175	omeg = tab_cntrl(tab0+6)
	176	g = tab_cntrl(tab0+7)
	177	mugaz = tab_cntrl(tab0+8)
	178	rcp = tab_cntrl(tab0+9)
	179	daysec = tab_cntrl(tab0+10)
	180	dtphys = tab_cntrl(tab0+11)
	181	c Informations about planet Mars for the physics only
	182	year_day = tab_cntrl(tab0+14)
	183	periheli = tab_cntrl(tab0+15)
	184	aphelie = tab_cntrl(tab0+16)
	185	peri_day = tab_cntrl(tab0+17)
	186	obliquit = tab_cntrl(tab0+18)
	187	c boundary layer and turbeulence
[224]	188	z0_default = tab_cntrl(tab0+19)
[38]	189	lmixmin = tab_cntrl(tab0+20)
	190	emin_turb = tab_cntrl(tab0+21)
	191	c optical properties of polar caps and ground emissivity
	192	albedice(1)= tab_cntrl(tab0+22)
	193	albedice(2)= tab_cntrl(tab0+23)
	194	emisice(1) = tab_cntrl(tab0+24)
	195	emisice(2) = tab_cntrl(tab0+25)
	196	emissiv = tab_cntrl(tab0+26)
	197	tauvis = tab_cntrl(tab0+27) ! dust opt. depth vis.
	198	iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
	199	iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
	200	dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north)
	201	dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south)
	202	c soil properties
	203	volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
	204	c-----------------------------------------------------------------------
	205	c Save some constants for later use (as routine arguments)
	206	c-----------------------------------------------------------------------
	207	p_omeg = omeg
	208	p_g = g
	209	p_mugaz = mugaz
	210	p_daysec = daysec
	211	p_rad=rad
	212
	213	ENDIF ! end of (nid = 0)
	214
	215	c-----------------------------------------------------------------------
	216	c Write physical constants to output before modifying them
	217	c-----------------------------------------------------------------------
	218
	219	6 FORMAT(a20,e15.6,e15.6)
	220	5 FORMAT(a20,f12.2,f12.2)
	221
	222	write(,) '*****************************************************'
	223	write(,) 'Reading tab_cntrl when calling tabfi before changes'
	224	write(,) '*****************************************************'
[1112]	225	write(*,5) '(1) = ngrid?',tab_cntrl(tab0+1)
[38]	226	write(*,5) '(2) lmax',tab_cntrl(tab0+2),real(lmax)
	227	write(*,5) '(3) day_ini',tab_cntrl(tab0+3),real(day_ini)
	228	write(*,5) '(5) rad',tab_cntrl(tab0+5),rad
	229	write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec
	230	write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg
	231	write(*,5) '(7) g',tab_cntrl(tab0+7),g
	232	write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz
	233	write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp
	234	write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys
	235
	236	write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day
	237	write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli
	238	write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie
	239	write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day
	240	write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit
	241
[224]	242	write(*,6) '(19) z0_default',tab_cntrl(tab0+19),z0_default
[38]	243	write(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb
	244	write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin
	245
	246	write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv
	247	write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1)
	248	write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2)
	249	write(*,5) '(22) albedice(1)',tab_cntrl(tab0+22),albedice(1)
	250	write(*,5) '(23) albedice(2)',tab_cntrl(tab0+23),albedice(2)
	251	write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
	252	write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
	253	write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
	254	write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
	255
	256	write(*,5) '(27) tauvis',tab_cntrl(tab0+27),tauvis
	257
	258	write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa
	259
	260	write(,)
	261	write(,) 'Lmodif in tabfi!!!!!!!',Lmodif
	262
	263	c-----------------------------------------------------------------------
	264	c Modifications...
[1130]	265	! NB: Modifying controls should only be done by newstart, and in seq mode
	266	if ((Lmodif.eq.1).and.is_parallel) then
	267	write(,) "tabfi: Error modifying tab_control should",
	268	& " only happen in serial mode (eg: by newstart)"
[2311]	269	call abort_physic(modname,
	270	& "tab_control modification not possible in parallel",1)
[1130]	271	endif
[38]	272	c-----------------------------------------------------------------------
	273
	274	IF(Lmodif.eq.1) then
	275
	276	write(,)
	277	write(,) 'Change values in tab_cntrl ? :'
	278	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	279	write(,) '(Current values given above)'
	280	write(,)
	281	write(,) '(3) day_ini : Initial day (=0 at Ls=0)'
[224]	282	write(,) '(19) z0 : default surface roughness (m)'
[38]	283	write(,) '(21) emin_turb : minimal energy (PBL)'
	284	write(,) '(20) lmixmin : mixing length (PBL)'
	285	write(,) '(26) emissiv : ground emissivity'
	286	write(,) '(24 et 25) emisice : CO2 ice max emissivity '
	287	write(,) '(22 et 23) albedice : CO2 ice cap albedos'
	288	write(,) '(31 et 32) iceradius : mean scat radius of CO2 snow'
	289	write(,) '(33 et 34) dtemisice : time scale for snow',
	290	& ' metamorphism'
	291	write(,) '(27) tauvis : mean dust vis. reference ',
	292	& 'opacity'
	293	write(,) '(35) volcapa : soil volumetric heat capacity'
	294	write(,) '(18) obliquit : planet obliquity (deg)'
	295	write(,) '(17) peri_day : perihelion date (sol since Ls=0)'
[669]	296	write(,) '( ) peri_ls : perihelion date (Ls since Ls=0)'
[38]	297	write(,) '(15) periheli : min. sun-mars dist (Mkm)'
	298	write(,) '(16) aphelie : max. sun-mars dist (Mkm)'
	299	write(,)
	300
	301
[552]	302	do ! neverending loop
[38]	303	write(,)
	304	write(,)
	305	write(,) 'Changes to perform ?'
	306	write(,) ' (enter keyword or return )'
	307	write(,)
	308	read(*,fmt='(a20)') modif
	309	if (modif(1:1) .eq. ' ') goto 999
	310
	311	write(,)
[552]	312	write(,) trim(modif) , ' : '
[38]	313
[552]	314	if (trim(modif) .eq. 'day_ini') then
[38]	315	write(,) 'current value:',day_ini
	316	write(,) 'enter new value:'
	317	101 read(,,iostat=ierr) day_ini
	318	if(ierr.ne.0) goto 101
	319	write(,) ' '
	320	write(,) 'day_ini (new value):',day_ini
	321
[552]	322	else if (trim(modif) .eq. 'z0') then
[224]	323	write(,) 'current value (m):',z0_default
	324	write(,) 'enter new value (m):'
	325	102 read(,,iostat=ierr) z0_default
[38]	326	if(ierr.ne.0) goto 102
	327	write(,) ' '
[224]	328	write(,) ' z0 (new value):',z0_default
[38]	329
[552]	330	else if (trim(modif) .eq. 'emin_turb') then
[38]	331	write(,) 'current value:',emin_turb
	332	write(,) 'enter new value:'
	333	103 read(,,iostat=ierr) emin_turb
	334	if(ierr.ne.0) goto 103
	335	write(,) ' '
	336	write(,) ' emin_turb (new value):',emin_turb
	337
[552]	338	else if (trim(modif) .eq. 'lmixmin') then
[38]	339	write(,) 'current value:',lmixmin
	340	write(,) 'enter new value:'
	341	104 read(,,iostat=ierr) lmixmin
	342	if(ierr.ne.0) goto 104
	343	write(,) ' '
	344	write(,) ' lmixmin (new value):',lmixmin
	345
[552]	346	else if (trim(modif) .eq. 'emissiv') then
[38]	347	write(,) 'current value:',emissiv
	348	write(,) 'enter new value:'
	349	105 read(,,iostat=ierr) emissiv
	350	if(ierr.ne.0) goto 105
	351	write(,) ' '
	352	write(,) ' emissiv (new value):',emissiv
	353
[552]	354	else if (trim(modif) .eq. 'emisice') then
[38]	355	write(,) 'current value emisice(1) North:',emisice(1)
	356	write(,) 'enter new value:'
	357	106 read(,,iostat=ierr) emisice(1)
	358	if(ierr.ne.0) goto 106
	359	write(,)
	360	write(,) ' emisice(1) (new value):',emisice(1)
	361	write(,)
	362
	363	write(,) 'current value emisice(2) South:',emisice(2)
	364	write(,) 'enter new value:'
	365	107 read(,,iostat=ierr) emisice(2)
	366	if(ierr.ne.0) goto 107
	367	write(,)
	368	write(,) ' emisice(2) (new value):',emisice(2)
	369
[552]	370	else if (trim(modif) .eq. 'albedice') then
[38]	371	write(,) 'current value albedice(1) North:',albedice(1)
	372	write(,) 'enter new value:'
	373	108 read(,,iostat=ierr) albedice(1)
	374	if(ierr.ne.0) goto 108
	375	write(,)
	376	write(,) ' albedice(1) (new value):',albedice(1)
	377	write(,)
	378
	379	write(,) 'current value albedice(2) South:',albedice(2)
	380	write(,) 'enter new value:'
	381	109 read(,,iostat=ierr) albedice(2)
	382	if(ierr.ne.0) goto 109
	383	write(,)
	384	write(,) ' albedice(2) (new value):',albedice(2)
	385
[552]	386	else if (trim(modif) .eq. 'iceradius') then
[38]	387	write(,) 'current value iceradius(1) North:',iceradius(1)
	388	write(,) 'enter new value:'
	389	110 read(,,iostat=ierr) iceradius(1)
	390	if(ierr.ne.0) goto 110
	391	write(,)
	392	write(,) ' iceradius(1) (new value):',iceradius(1)
	393	write(,)
	394
	395	write(,) 'current value iceradius(2) South:',iceradius(2)
	396	write(,) 'enter new value:'
	397	111 read(,,iostat=ierr) iceradius(2)
	398	if(ierr.ne.0) goto 111
	399	write(,)
	400	write(,) ' iceradius(2) (new value):',iceradius(2)
	401
[552]	402	else if (trim(modif) .eq. 'dtemisice') then
[38]	403	write(,) 'current value dtemisice(1) North:',dtemisice(1)
	404	write(,) 'enter new value:'
	405	112 read(,,iostat=ierr) dtemisice(1)
	406	if(ierr.ne.0) goto 112
	407	write(,)
	408	write(,) ' dtemisice(1) (new value):',dtemisice(1)
	409	write(,)
	410
	411	write(,) 'current value dtemisice(2) South:',dtemisice(2)
	412	write(,) 'enter new value:'
	413	113 read(,,iostat=ierr) dtemisice(2)
	414	if(ierr.ne.0) goto 113
	415	write(,)
	416	write(,) ' dtemisice(2) (new value):',dtemisice(2)
	417
[552]	418	else if (trim(modif) .eq. 'tauvis') then
[38]	419	write(,) 'current value:',tauvis
	420	write(,) 'enter new value:'
	421	114 read(,,iostat=ierr) tauvis
	422	if(ierr.ne.0) goto 114
	423	write(,)
	424	write(,) ' tauvis (new value):',tauvis
	425
[552]	426	else if (trim(modif) .eq. 'obliquit') then
[38]	427	write(,) 'current value:',obliquit
	428	write(,) 'obliquit should be 25.19 on current Mars'
	429	write(,) 'enter new value:'
	430	115 read(,,iostat=ierr) obliquit
	431	if(ierr.ne.0) goto 115
	432	write(,)
	433	write(,) ' obliquit (new value):',obliquit
	434
[552]	435	else if (trim(modif) .eq. 'peri_day') then
[38]	436	write(,) 'current value:',peri_day
[669]	437	write(,) 'peri_day should be 485 sols on current Mars'
[38]	438	write(,) 'enter new value:'
	439	116 read(,,iostat=ierr) peri_day
	440	if(ierr.ne.0) goto 116
	441	write(,)
	442	write(,) ' peri_day (new value):',peri_day
[669]	443
	444	else if (trim(modif) .eq. 'peri_ls') then
	445	write(,) 'peri_ls value is not stored in start files,'
	446	write(,) 'but it should be 251 degrees on current Mars'
	447	write(,) '(peri_day should be 485 sols on current Mars)'
	448	write(,) 'enter new value:'
	449	1160 read(,,iostat=ierr) peri_ls
	450	if(ierr.ne.0) goto 1160
	451	write(,)
[672]	452	write(,) 'peri_ls asked:',peri_ls
	453	write(,) 'for aphelion =',aphelie
	454	write(,) 'perihelion =',periheli
	455	write(,) 'and',year_day,'sols/year'
	456	call lsp2solp(peri_ls,peri_day,aphelie,periheli,year_day)
	457	write(,) 'peri_day (new value):',peri_day
[38]	458
[669]	459
[552]	460	else if (trim(modif) .eq. 'periheli') then
[38]	461	write(,) 'current value:',periheli
	462	write(,) 'perihelion should be 206.66 on current Mars'
	463	write(,) 'enter new value:'
	464	117 read(,,iostat=ierr) periheli
	465	if(ierr.ne.0) goto 117
	466	write(,)
	467	write(,) ' periheli (new value):',periheli
	468
[552]	469	else if (trim(modif) .eq. 'aphelie') then
[38]	470	write(,) 'current value:',aphelie
	471	write(,) 'aphelion should be 249.22 on current Mars'
	472	write(,) 'enter new value:'
	473	118 read(,,iostat=ierr) aphelie
	474	if(ierr.ne.0) goto 118
	475	write(,)
	476	write(,) ' aphelie (new value):',aphelie
	477
[552]	478	else if (trim(modif) .eq. 'volcapa') then
[38]	479	write(,) 'current value:',volcapa
	480	write(,) 'enter new value:'
	481	119 read(,,iostat=ierr) volcapa
	482	if(ierr.ne.0) goto 119
	483	write(,)
	484	write(,) ' volcapa (new value):',volcapa
	485
	486	endif
[552]	487	enddo ! of do ! neverending loop
[38]	488
	489	999 continue
	490
	491	c-----------------------------------------------------------------------
	492	c Write values of physical constants after modifications
	493	c-----------------------------------------------------------------------
	494
	495	write(,) '*****************************************************'
	496	write(,) 'Reading tab_cntrl when calling tabfi AFTER changes'
	497	write(,) '*****************************************************'
[1112]	498	write(*,5) '(1) = ngrid?',tab_cntrl(tab0+1)
[38]	499	write(*,5) '(2) lmax',tab_cntrl(tab0+2),real(lmax)
	500	write(*,5) '(3) day_ini',tab_cntrl(tab0+3),real(day_ini)
	501	write(*,5) '(5) rad',tab_cntrl(tab0+5),rad
	502	write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec
	503	write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg
	504	write(*,5) '(7) g',tab_cntrl(tab0+7),g
	505	write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz
	506	write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp
	507	write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys
	508
	509	write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day
	510	write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli
	511	write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie
	512	write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day
	513	write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit
	514
[224]	515	write(*,6) '(19) z0_default',tab_cntrl(tab0+19),z0_default
[38]	516	write(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb
	517	write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin
	518
	519	write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv
	520	write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1)
	521	write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2)
	522	write(*,5) '(22) albedice(1)',tab_cntrl(tab0+22),albedice(1)
	523	write(*,5) '(23) albedice(2)',tab_cntrl(tab0+23),albedice(2)
	524	write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
	525	write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
	526	write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
	527	write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
	528
	529	write(*,5) '(27) tauvis',tab_cntrl(tab0+27),tauvis
	530
	531	write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa
	532
	533	write(,)
	534	write(,)
	535
[1524]	536	ENDIF ! of if (Lmodif == 1)
[38]	537
	538	c-----------------------------------------------------------------------
	539	c Case when using a start file from before March 1996 (without iceradius...
	540	c-----------------------------------------------------------------------
	541	if (iceradius(1).eq.0) then
	542	iceradius(1) = 100.e-6
	543	iceradius(2) = 100.e-6
	544	dtemisice(1) = 0.4
	545	dtemisice(2) = 0.4
	546	write (,) ' tabfi: WARNING : old initialisation file'
	547	write (,) 'iceradius set to',iceradius(1),iceradius(2)
	548	write (,) 'dtemisice set to',dtemisice(1),dtemisice(2)
	549	end if
	550
	551	c-----------------------------------------------------------------------
	552	end
[669]	553
	554
	555
	556
	557	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	558	! gives sol at perihelion for ls at perihelion (for precession cycles)
[672]	559	subroutine lsp2solp(lsp,solp,aphelie,periheli,year_day)
[669]	560
	561	implicit none
	562	! Arguments:
	563	real lsp ! Input: ls at perihelion
	564	real solp ! Output: sol at perihelion
[672]	565	real aphelie,periheli,year_day ! Input: parameters
	566
[669]	567	! Local:
	568	double precision zx0 ! eccentric anomaly at Ls=0
	569	double precision e_elips
	570	double precision pi,degrad
	571
	572	parameter (pi=3.14159265358979d0)
	573	parameter (degrad=57.2957795130823d0)
[672]	574
	575	e_elips=(aphelie-periheli)/(aphelie+periheli)
[669]	576	zx0 = -2.0datan(dtan(0.5lsp/degrad)
	577	. *dsqrt((1.-e_elips)/(1.+e_elips)))
	578	if (zx0 .le. 0.) zx0 = zx0 + 2.*pi
	579
	580	solp = year_day(1.-(zx0-e_elipsdsin(zx0))/(2.*pi))
	581
	582
	583	end subroutine lsp2solp
	584
	585
	586

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

Download in other formats:

Original Format