Context Navigation

tabfi.F @ 1655

Last change on this file since 1655 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

File size: 22.8 KB

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

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

Download in other formats:

Original Format