Context Navigation

tabfi.F @ 1431

Last change on this file since 1431 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format