Context Navigation

tabfi.F @ 1246

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

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

Download in other formats:

Original Format