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

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
53	use planete_h
54	implicit none
55
56	#include "netcdf.inc"
57
58	c-----------------------------------------------------------------------
59	c Declarations
60	c-----------------------------------------------------------------------
61
62	c Arguments
63	c ---------
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
69
70	c Variables
71	c ---------
72	INTEGER :: nvarid
73	REAL :: peri_ls
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
80	LOGICAL :: found
81
82	write(,)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
83
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----------------------------
112	z0_default = 1.e-2 ! surface roughness (m) ~0.01
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
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
162
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
170	c
171	c Initialization of some physical constants
172	c informations on physics grid
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
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
196	z0_default = tab_cntrl(tab0+19)
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(,) '*****************************************************'
233	write(*,5) '(1) = ngrid?',tab_cntrl(tab0+1)
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
250	write(*,6) '(19) z0_default',tab_cntrl(tab0+19),z0_default
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...
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
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)'
289	write(,) '(19) z0 : default surface roughness (m)'
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)'
303	write(,) '( ) peri_ls : perihelion date (Ls since Ls=0)'
304	write(,) '(15) periheli : min. sun-mars dist (Mkm)'
305	write(,) '(16) aphelie : max. sun-mars dist (Mkm)'
306	write(,)
307
308
309	do ! neverending loop
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(,)
319	write(,) trim(modif) , ' : '
320
321	if (trim(modif) .eq. 'day_ini') then
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
329	else if (trim(modif) .eq. 'z0') then
330	write(,) 'current value (m):',z0_default
331	write(,) 'enter new value (m):'
332	102 read(,,iostat=ierr) z0_default
333	if(ierr.ne.0) goto 102
334	write(,) ' '
335	write(,) ' z0 (new value):',z0_default
336
337	else if (trim(modif) .eq. 'emin_turb') then
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
345	else if (trim(modif) .eq. 'lmixmin') then
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
353	else if (trim(modif) .eq. 'emissiv') then
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
361	else if (trim(modif) .eq. 'emisice') then
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
377	else if (trim(modif) .eq. 'albedice') then
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
393	else if (trim(modif) .eq. 'iceradius') then
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
409	else if (trim(modif) .eq. 'dtemisice') then
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
425	else if (trim(modif) .eq. 'tauvis') then
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
433	else if (trim(modif) .eq. 'obliquit') then
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
442	else if (trim(modif) .eq. 'peri_day') then
443	write(,) 'current value:',peri_day
444	write(,) 'peri_day should be 485 sols on current Mars'
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
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(,)
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
465
466
467	else if (trim(modif) .eq. 'periheli') then
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
476	else if (trim(modif) .eq. 'aphelie') then
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
485	else if (trim(modif) .eq. 'volcapa') then
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
494	enddo ! of do ! neverending loop
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(,) '*****************************************************'
505	write(*,5) '(1) = ngrid?',tab_cntrl(tab0+1)
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
522	write(*,6) '(19) z0_default',tab_cntrl(tab0+19),z0_default
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
560
561
562
563
564	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
565	! gives sol at perihelion for ls at perihelion (for precession cycles)
566	subroutine lsp2solp(lsp,solp,aphelie,periheli,year_day)
567
568	implicit none
569	! Arguments:
570	real lsp ! Input: ls at perihelion
571	real solp ! Output: sol at perihelion
572	real aphelie,periheli,year_day ! Input: parameters
573
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)
581
582	e_elips=(aphelie-periheli)/(aphelie+periheli)
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