Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

tabfi.F @ 3572

Last change on this file since 3572 was 3175, checked in by emillour, 18 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 18.2 KB

Rev	Line
[3175]	1	c=======================================================================
	2	SUBROUTINE tabfi(nid,Lmodif,tab0,day_ini,lmax,p_rad,
	3	. p_omeg,p_g,p_cpp,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	! to use 'getin'
	46	use ioipsl_getincom
	47	use planet_h
	48	use comgeomfi_h
	49	use comsoil_h
	50	implicit none
	51
	52	#include "dimensions.h"
	53	#include "dimphys.h"
	54	#include "comcstfi.h"
	55	#include "surfdat.h"
	56	#include "netcdf.inc"
	57	#include "callkeys.h"
	58
	59	c-----------------------------------------------------------------------
	60	c Declarations
	61	c-----------------------------------------------------------------------
	62
	63	c Arguments
	64	c ---------
	65	INTEGER nid,nvarid,tab0
	66	INTEGER*4 day_ini
	67	INTEGER Lmodif
	68	INTEGER lmax
	69	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
	70	REAL rm !TB
	71
	72	c Variables
	73	c ---------
	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 modif20,choice1
	80	REAL eccpal,sma ! paleo
	81	c-----------------------------------------------------------------------
	82	c Initialization of various physical constants to defaut values (nid = 0 case)
	83	c-----------------------------------------------------------------------
	84	c TB : modif : suppression de la boucle nid
	85	c IF (nid.eq.0) then
	86
	87	c Reference pressure
	88	c-------------------------------------
	89	c pressrf = 1. ! Pression de reference (Pa) ~1
	90
	91	c Default (Pluton) parameters for the dynamics and physics
	92	c----------------------------------------------------------
	93	rad=1187000 !old= 1173000. ! radius of Pluton (m)
	94	daysec=551854.08 ! length of a sol (s) ~551837 s
	95	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
	96	g=0.617189 ! gravity (m.s-2) GM/r2 old: 0.65
	97	mugaz=28. ! Molar mass of the atmosphere (g.mol-1) ~28
	98	rcp=.2853 ! = r/cp ~0.2853
	99
	100	c Default (Pluton) parameters for physics only
	101	c----------------------------------------------
	102	year_day = 14178.30343 ! length of year (sols) ~14182.245
	103	periheli = 4436.82 ! min. Star-Planet distance (Mkm) ~4436
	104	aphelie = 7375.93 ! max. Star-Planet distance (Mkm) ~7375
	105	peri_day = 87.362 ! date of perihelion (sols since N. spring)
	106	obliquit=119.591 ! Obliquity of the planet (deg) old ~119.6
	107	tpal=0. ! time (Myr) to start paleoclimate run (e.g -10)
	108	adjust=0. ! N2 ice albedo adjustment
	109
	110	c Boundary layer and turbulence
	111	c----------------------------
	112	z0 = 1.e-2 ! surface roughness (m) ~0.01
	113
	114	c soil properties
	115	volcapa = 1.e6 ! soil volumetric heat capacity (in comsoil.h)
	116
	117	c-----------------------------------------------------------------------
	118	c Initialization of physical constants by reading array tab_cntrl(:)
	119	c which contains these parameters (nid != 0 case)
	120	c-----------------------------------------------------------------------
	121	c Read 'controle' array
	122	c
	123	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	124	IF (ierr .NE. NF_NOERR) THEN
	125	PRINT*, "Tabfi: Could not find <controle> data"
	126	CALL abort
	127	ENDIF
	128	#ifdef NC_DOUBLE
	129	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	130	#else
	131	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	132	#endif
	133	IF (ierr .NE. NF_NOERR) THEN
	134	PRINT*, "Tabfi: Failed reading <controle> array"
	135	CALL abort
	136	ENDIF
	137
	138	c Write physical constants to output before modifying them
	139	c-----------------------------------------------------------------------
	140
	141	6 FORMAT(a20,e15.6,e15.6)
	142	5 FORMAT(a20,f12.2,f12.2)
	143
	144	write(,) '*****************************************************'
	145	write(,) 'Reading tab_cntrl when calling tabfi before changes'
	146	write(,) '*****************************************************'
	147	write(,) ' VALUE IN NC DEFAULT VALUE'
	148	write(*,5) '(1) = ngridmx?',tab_cntrl(tab0+1),float(ngridmx)
	149	write(*,5) '(2) lmax',tab_cntrl(tab0+2),float(lmax)
	150	write(*,5) '(3) day_ini',tab_cntrl(tab0+3),float(day_ini)
	151	write(*,5) '(5) rad',tab_cntrl(tab0+5),rad
	152	write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec
	153	write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg
	154	write(*,5) '(7) g',tab_cntrl(tab0+7),g
	155	write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz
	156	write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp
	157	write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys
	158
	159	write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day
	160	write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli
	161	write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie
	162	write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day
	163	write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit
	164
	165	write(*,6) '(19) z0',tab_cntrl(tab0+19),z0
	166	write(*,5) '(20) tpal',tab_cntrl(tab0+20),tpal
	167	write(*,5) '(21) adjust',tab_cntrl(tab0+21),adjust
	168
	169	write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa
	170
	171	write(,)
	172
	173	c
	174	c Initialization of some physical constants
	175	c informations on physics grid
	176	if(ngridmx.ne.tab_cntrl(tab0+1)) then
	177	print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngridmx'
	178	print*,tab_cntrl(tab0+1),ngridmx
	179	endif
	180	lmax = nint(tab_cntrl(tab0+2))
	181	day_ini = tab_cntrl(tab0+3)
	182	time = tab_cntrl(tab0+4)
	183	c Informations about planet for dynamics and physics
	184	rad = tab_cntrl(tab0+5)
	185	omeg = tab_cntrl(tab0+6)
	186	g = tab_cntrl(tab0+7)
	187	mugaz = tab_cntrl(tab0+8)
	188	rcp = tab_cntrl(tab0+9)
	189	daysec = tab_cntrl(tab0+10)
	190	dtphys = tab_cntrl(tab0+11)
	191	c Informations about planet for the physics only
	192	year_day = tab_cntrl(tab0+14)
	193	periheli = tab_cntrl(tab0+15)
	194	aphelie = tab_cntrl(tab0+16)
	195	peri_day = tab_cntrl(tab0+17)
	196	obliquit = tab_cntrl(tab0+18)
	197	c boundary layer and turbeulence
	198	z0 = tab_cntrl(tab0+19)
	199	c for paleoclimate
	200	tpal = tab_cntrl(tab0+20)
	201	adjust = tab_cntrl(tab0+21) ! for Triton albedo adjustment
	202	c soil properties
	203	volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
	204
	205	c-----------------------------------------------------------------------
	206	c Save some constants for later use (as routine arguments)
	207	c-----------------------------------------------------------------------
	208	p_omeg = omeg
	209	p_g = g
	210	rm= 8.314511E+0 *1000.E+0/mugaz
	211	p_cpp = rm/rcp
	212	p_mugaz = mugaz
	213	p_daysec = daysec
	214	p_rad=rad
	215	write(,),'In particular in the nc file: rad, g = ',p_rad,p_g
	216	write(,),'cpp calculated from rcp is : ',p_cpp
	217
	218	c ENDIF ! end of (nid = 0)
	219
	220	c-----------------------------------------------------------------------
	221	c Modifications...
	222	c-----------------------------------------------------------------------
	223	IF(Lmodif.eq.1) then
	224
	225	write(,)
	226	write(,) 'Change values in tab_cntrl ? :'
	227	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	228	write(,) '(Current values given above)'
	229	write(,)
	230	write(,) '(2) lmax'
	231	write(,) '(3) day_ini : Initial day (=0 at Ls=0)'
	232	write(,) '(5) rad : radius of the planet (m)'
	233	write(,) '(6) omeg : planet rotation rate (rad/s)'
	234	write(,) '(7) g : gravity (m/s2)'
	235	write(,) '(8) mugaz : molecular mass of atm (g/mol)'
	236	write(,) '(9) rcp : r/Cp'
	237	write(,) '(10) daysec : length of a sol (s)'
	238	write(,) '(14) year_day : length of year (in sols)'
	239	write(,) '(15) periheli : min. sun-pluto dist (Mkm)'
	240	write(,) '(16) aphelie : max. sun-pluto dist (Mkm)'
	241	write(,) '(17) peri_day : perihelion date (sol since Ls=0)'
	242	write(,) '(18) obliquit : planet obliquity (deg)'
	243	write(,) '(19) z0 : surface roughness (m)'
	244	write(,) '(20) tpal: paleoclimatic date in million years'
	245	write(,) '(21) adjust: N2 ice albedo adjustment'
	246	write(,) '(35) volcapa : soil volumetric heat capacity'
	247
	248	write(,)
	249
	250	do while(modif(1:1).ne.'hello')
	251	write(,)
	252	write(,)
	253	write(,) 'Changes to perform ?'
	254	write(,) ' (enter keyword or return )'
	255	write(,)
	256	read(*,fmt='(a20)') modif
	257	if (modif(1:1) .eq. ' ') goto 999
	258
	259	write(,)
	260	write(,) modif(1:len_trim(modif)) , ' : '
	261
	262	if (modif(1:len_trim(modif)) .eq. 'day_ini') then
	263	write(,) 'current value:',day_ini
	264	write(,) 'enter new value:'
	265	101 read(,,iostat=ierr) day_ini
	266	if(ierr.ne.0) goto 101
	267	write(,) ' '
	268	write(,) 'day_ini (new value):',day_ini
	269
	270	else if (modif(1:len_trim(modif)) .eq. 'z0') then
	271	write(,) 'current value:',z0
	272	write(,) 'enter new value:'
	273	102 read(,,iostat=ierr) z0
	274	if(ierr.ne.0) goto 102
	275	write(,) ' '
	276	write(,) ' z0 (new value):',z0
	277
	278	else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
	279	write(,) 'current value:',obliquit
	280	write(,) 'obliquit should be ... on current Pluto'
	281	write(,) 'enter new value:'
	282	115 read(,,iostat=ierr) obliquit
	283	if(ierr.ne.0) goto 115
	284	write(,)
	285	write(,) ' obliquit (new value):',obliquit
	286
	287	else if (modif(1:len_trim(modif)) .eq. 'peri_day') then
	288	write(,) 'current value:',peri_day
	289	write(,) 'peri_day should be 72 on current Pluto'
	290	write(,) 'enter new value:'
	291	116 read(,,iostat=ierr) peri_day
	292	if(ierr.ne.0) goto 116
	293	write(,)
	294	write(,) ' peri_day (new value):',peri_day
	295
	296	else if (modif(1:len_trim(modif)) .eq. 'periheli') then
	297	write(,) 'current value:',periheli,' e6 km'
	298	write(,) 'perihelion should be 4436 on current Pluto'
	299	write(,) 'enter new value:'
	300	117 read(,,iostat=ierr) periheli
	301	if(ierr.ne.0) goto 117
	302	write(,)
	303	write(,) ' periheli (new value):',periheli
	304
	305	else if (modif(1:len_trim(modif)) .eq. 'aphelie') then
	306	write(,) 'current value:',aphelie,' e6 km'
	307	write(,) 'aphelion should be 7375 on current Pluto'
	308	write(,) 'enter new value:'
	309	118 read(,,iostat=ierr) aphelie
	310	if(ierr.ne.0) goto 118
	311	write(,)
	312	write(,) ' aphelie (new value):',aphelie
	313
	314	else if (modif(1:len_trim(modif)) .eq. 'volcapa') then
	315	write(,) 'current value:',volcapa
	316	write(,) 'enter new value:'
	317	119 read(,,iostat=ierr) volcapa
	318	if(ierr.ne.0) goto 119
	319	write(,)
	320	write(,) ' volcapa (new value):',volcapa
	321
	322	else if (modif(1:len_trim(modif)) .eq. 'year_day') then
	323	write(,) 'current value:',year_day
	324	write(,) 'enter new value:'
	325	120 read(,,iostat=ierr) year_day
	326	if(ierr.ne.0) goto 120
	327	write(,)
	328	write(,) ' year_day (new value):',year_day
	329
	330	else if (modif(1:len_trim(modif)).eq.'rad') then
	331	write(,) 'current value:',rad
	332	write(,) 'enter new value:'
	333	121 read(,,iostat=ierr) rad
	334	if(ierr.ne.0) goto 121
	335	write(,)
	336	write(,) ' rad (new value):',rad
	337
	338	else if (modif(1:len_trim(modif)).eq.'omeg') then
	339	write(,) 'current value:',omeg
	340	write(,) 'enter new value:'
	341	122 read(,,iostat=ierr) omeg
	342	if(ierr.ne.0) goto 122
	343	write(,)
	344	write(,) ' omeg (new value):',omeg
	345
	346	else if (modif(1:len_trim(modif)).eq.'g') then
	347	write(,) 'current value:',g
	348	write(,) 'enter new value:'
	349	123 read(,,iostat=ierr) g
	350	if(ierr.ne.0) goto 123
	351	write(,)
	352	write(,) ' g (new value):',g
	353
	354	else if (modif(1:len_trim(modif)).eq.'mugaz') then
	355	write(,) 'current value:',mugaz
	356	write(,) 'enter new value:'
	357	124 read(,,iostat=ierr) mugaz
	358	if(ierr.ne.0) goto 124
	359	write(,)
	360	write(,) ' mugaz (new value):',mugaz
	361	r=8.314/(mugaz/1000.0)
	362	rm=r
	363	write(,) ' r (new value):',r
	364
	365	else if (modif(1:len_trim(modif)).eq.'rcp') then
	366	write(,) 'current value:',rcp
	367	write(,) 'enter new value:'
	368	125 read(,,iostat=ierr) rcp
	369	if(ierr.ne.0) goto 125
	370	write(,)
	371	write(,) ' rcp (new value):',rcp
	372	cpp=rm/rcp
	373	write(,) ' cpp (new value):',cpp
	374
	375	else if (modif(1:len_trim(modif)).eq.'daysec') then
	376	write(,) 'current value:',daysec
	377	write(,) 'enter new value:'
	378	126 read(,,iostat=ierr) daysec
	379	if(ierr.ne.0) goto 126
	380	write(,)
	381	write(,) ' daysec (new value):',daysec
	382	write(,) ' want to update omeg (y/n)?'
	383	128 read(,,iostat=ierr) choice
	384	if(ierr.ne.0) goto 128
	385	if (choice.eq.'y') then
	386	omeg=2.*3.14159265/daysec
	387	write(,) ' omeg (new value):',omeg
	388	endif
	389
	390	else if (modif(1:len_trim(modif)).eq.'lmax') then
	391	write(,) 'current value:',lmax
	392	write(,) 'enter new value:'
	393	127 read(,,iostat=ierr) lmax
	394	if(ierr.ne.0) goto 127
	395	write(,)
	396	write(,) ' lmax (new value):',lmax
	397
	398	else if (modif(1:len_trim(modif)).eq.'tpal') then
	399	write(,) 'current value:',tpal
	400	write(,) 'enter new value:'
	401	129 read(,,iostat=ierr) tpal
	402	if(ierr.ne.0) goto 129
	403	write(,)
	404	write(,) ' tpal (new value):',tpal
	405	write(,) 'Modify obli / peri_day / ecc accordingly (y/n)?'
	406	130 read(,,iostat=ierr) choice
	407	if(ierr.ne.0) goto 130
	408	if (choice.eq.'y') then
	409	c Initialising orbital parameters:
	410	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	411	c Updating obliquity and orbital parameters:
	412	call calcmilank(tpal,obliquit,peri_day,eccpal)
	413	sma=(periheli+aphelie)/2.
	414	periheli=sma*(1-eccpal)
	415	aphelie=sma*(1+eccpal)
	416	endif
	417
	418	else if (modif(1:len_trim(modif)).eq.'adjust') then
	419	write(,) 'current value:',adjust
	420	write(,) 'enter new value:'
	421	131 read(,,iostat=ierr) adjust
	422	if(ierr.ne.0) goto 131
	423	write(,)
	424	write(,) ' adjust (new value):',adjust
	425
	426
	427	endif
	428	enddo ! of do while(modif(1:1).ne.'hello')
	429
	430	999 continue
	431
	432	c-----------------------------------------------------------------------
	433	c Write values of physical constants after modifications
	434	c-----------------------------------------------------------------------
	435
	436	write(,) '*****************************************************'
	437	write(,) 'Reading tab_cntrl when calling tabfi AFTER changes'
	438	write(,) '*****************************************************'
	439	write(,) ' OLD VALUE IN NC NEW VALUE IN NC'
	440	write(*,5) '(1) = ngridmx?',tab_cntrl(tab0+1),float(ngridmx)
	441	write(*,5) '(2) lmax',tab_cntrl(tab0+2),float(lmax)
	442	write(*,5) '(3) day_ini',tab_cntrl(tab0+3),float(day_ini)
	443	write(*,5) '(5) rad',tab_cntrl(tab0+5),rad
	444	write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec
	445	write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg
	446	write(*,5) '(7) g',tab_cntrl(tab0+7),g
	447	write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz
	448	write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp
	449	write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys
	450
	451	write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day
	452	write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli
	453	write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie
	454	write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day
	455	write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit
	456
	457	write(*,6) '(19) z0',tab_cntrl(tab0+19),z0
	458	write(*,5) '(20) tpal',tab_cntrl(tab0+20),tpal
	459	write(*,5) '(21) adjust',tab_cntrl(tab0+21),adjust
	460
	461	write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa
	462
	463	write(,)
	464	write(,)
	465
	466	ENDIF ! of if (Lmodif == 1)
	467
	468	c-----------------------------------------------------------------------
	469	c Save some constants for later use (as routine arguments)
	470	c-----------------------------------------------------------------------
	471	p_omeg = omeg
	472	p_g = g
	473	rm= 8.314511E+0 *1000.E+0/mugaz
	474	p_cpp = rm/rcp
	475	p_mugaz = mugaz
	476	p_daysec = daysec
	477	p_rad=rad
	478	write(,),'NOW in the nc file: rad, g, omeg = ',p_rad,p_g,p_omeg
	479	write(,),'r is = ',rm,' mugaz is ',mugaz
	480	write(,),'cpp calculated from r and rcp is : ',p_cpp
	481	write(,),'daysec is : ',daysec
	482
	483	end

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

Download in other formats:

Original Format