Context Navigation

rcm1d.F @ 888

Last change on this file since 888 was 887, checked in by slebonnois, 13 years ago
SL: add rcm1d tool in Venus and Titan physics to run the code in 1D, and associated small modifications
File size: 13.4 KB

Line
1
2	PROGRAM rcm1d
3
4	USE infotrac
5	use control_mod
6	use comgeomphy
7	IMPLICIT NONE
8
9	c=======================================================================
10	c subject:
11	c --------
12	c PROGRAM useful to run physical part of the venusian GCM in a 1D column
13	c
14	c Can be compiled with a command like (e.g. for 50 layers)
15	c "makelmdz -p venus -d 50 rcm1d"
16
17	c It requires the files "rcm1d.def" "physiq.def"
18	c and a file describing the sigma layers (e.g. "z2sig.def")
19	c
20	c author: Frederic Hourdin, R.Fournier,F.Forget (original Mars version)
21	c ------- Sebastien Lebonnois (Venus version)
22	c
23	c=======================================================================
24
25	#include "dimensions.h"
26	#include "dimsoil.h"
27	#include "comcstfi.h"
28	#include "comvert.h"
29	#include "netcdf.inc"
30	#include "logic.h"
31	#include "clesphys.h"
32	#include "iniprint.h"
33	#include "tabcontrol.h"
34
35	c --------------------------------------------------------------
36	c Declarations
37	c --------------------------------------------------------------
38	c
39	INTEGER unit ! unite de lecture de "rcm1d.def"
40	INTEGER unitstart ! unite d'ecriture de "startphy.nc"
41	INTEGER nlayer,nlevel,nsoil,ndt
42	INTEGER ilayer,ilevel,isoil,idt,iq
43	LOGICAl firstcall,lastcall
44	c
45	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
46	REAL day ! date durant le run
47	REAL time ! time (0<time<1 ; time=0.5 a midi)
48	REAL play(llm) ! Pressure at the middle of the layers (Pa)
49	REAL plev(llm+1) ! intermediate pressure levels (pa)
50	REAL psurf,tsurf(1)
51	REAL u(llm),v(llm) ! zonal, meridional wind
52	REAL gru,grv ! prescribed "geostrophic" background wind
53	REAL temp(llm) ! temperature at the middle of the layers
54	REAL q(llm,nqtot) ! tracer mixing ratio (e.g. kg/kg)
55	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
56	REAL zlay(llm) ! altitude estimee dans les couches (km)
57	REAL long(1),lati(1),area(1)
58	REAL cufi(1),cvfi(1)
59	REAL phisfi(1),albedo(1)
60	REAL solsw(1),sollwdown(1),dlw(1),radsol(1)
61	REAL zmea(1), zstd(1)
62	REAL zsig(1), zgam(1), zthe(1)
63	REAL zpic(1), zval(1)
64
65	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
66	REAL du(llm),dv(llm),dtemp(llm)
67	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
68	REAL dpsurf
69	REAL dq(llm,nqtot)
70	REAL dqdyn(llm,nqtot)
71
72	c Various intermediate variables
73	REAL zls
74	REAL phi(llm),s(llm)
75	REAL pk(llm),pks, w(llm)
76	INTEGER l, ierr, aslun
77	REAL tmp1(0:llm),tmp2(0:llm)
78
79	character*2 str2
80
81	c normalement dans dyn3d/comconst.h
82	COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
83	REAL cppdyn,nu_venus,t0_venus
84
85	c=======================================================================
86
87	c=======================================================================
88	c INITIALISATION
89	c=======================================================================
90
91	lunout = 6
92
93	c ------------------------------------------------------
94	c Constantes prescrites ICI
95	c ------------------------------------------------------
96
97	pi=2.E+0*asin(1.E+0)
98
99	c Constante de la Planete Venus
100	c -----------------------------
101	planet_type = "venus"
102	rad=6051300. ! rayon de Venus (m) ~6051300 m
103	daysec= 1.0087e7 ! duree du sol (s) ~1.e7 s
104	omeg=4.*asin(1.)/19.4141e6 ! vitesse de rotation (rad.s-1)
105	g= 8.87 ! gravite (m.s-2) ~8.87
106	mugaz=43.44 ! Masse molaire de l'atm (g.mol-1) ~43.44
107	! ADAPTATION GCM POUR CP(T)
108	! VENUS: Cp(T) = cpp*(T/T0)^nu
109	! avec T0=460. et nu=0.35
110	cpp=1.0e3
111	! cpp=9.0e2 ! version constante
112	cppdyn=cpp
113	r= 8.314511E+0 *1000.E+0/mugaz
114	rcp= r/cpp
115
116	c ------------------------------------------------------
117	c Lecture des parametres dans "rcm1d.def"
118	c ------------------------------------------------------
119
120	c Opening parameters file "rcm1d.def"
121	c ---------------------------------------
122	unit =97
123	OPEN(unit,file='rcm1d.def',status='old',form='formatted'
124	. ,iostat=ierr)
125
126	IF(ierr.ne.0) THEN
127	write(,) 'Problem to open "rcm1d.def'
128	write(,) 'Is it there ?'
129	stop
130	END IF
131
132	c Date et heure locale du debut du run
133	c ------------------------------------
134	c Date (en sols depuis le solstice de printemps) du debut du run
135	day0 = 0
136	PRINT *,'date de depart ?'
137	READ(unit,*) day0
138	day=REAL(day0)
139	PRINT *,day0
140	c Heure de demarrage
141	PRINT *,'heure de debut de simulation (entre 0 et 24) ?'
142	READ(unit,*) time
143	time=time/24.E+0
144
145	c Discretisation (Definition de la grille et des pas de temps)
146	c --------------
147	c
148	nlayer=llm
149	nlevel=nlayer+1
150	nsoil=nsoilmx
151	PRINT *,'nombre de pas de temps par jour ?'
152	READ(unit,*) day_step
153	print*,day_step
154
155	c PRINT *,'nombre d appel au rayonnement par jour ?'
156	c READ(unit,*) nbapp_rad
157	c print*,nbapp_rad
158	c LU DANS PHYSIQ.DEF...
159	nbapp_rad = 1000.
160
161	PRINT *,'nombre de jours simules ?'
162	READ(unit,*) ndt
163	print*,ndt
164
165	ndt=ndt*day_step
166	dtphys=daysec/day_step
167
168	c Pression de surface sur la planete
169	c ------------------------------------
170	c
171	PRINT *,'pression au sol'
172	READ(unit,*) psurf
173	PRINT *,psurf
174	c Pression de reference ! voir dyn3d/etat0_venus
175	c pa = 5.e4
176	pa = 1.e6
177	preff = 9.2e6 ! 92 bars
178	c preff = psurf
179
180	c latitude/longitude
181	c -------------------
182	PRINT *,'latitude en degres ?'
183	READ(unit,*) lati(1)
184	PRINT *,lati(1)
185	long(1)=0.E+0
186
187	c Initialisation albedo
188	c ----------------------
189	c ne sert pas ici...
190	albedo(1)=0.1
191	c PRINT *,'Albedo du sol nu ?'
192	c READ(unit,*) albedo(1)
193	c PRINT *,albedo(1)
194
195	c Initialisation speciales "physiq"
196	c ---------------------------------
197
198	CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
199	call initcomgeomphy
200	call infotrac_init
201	call ini_cpdet
202
203	c la surface de chaque maille est inutile en 1D --->
204	area(1)=1.E+0
205	c de meme ?
206	cufi(1)=1.E+0
207	cvfi(1)=1.E+0
208
209	c le geopotentiel au sol est inutile en 1D car tout est controle
210	c par la pression de surface --->
211	phisfi(1)=0.E+0
212
213	CALL iniphysiq(1,llm,daysec,day0,dtphys,
214	. lati,long,area,cufi,cvfi,rad,g,r,cpp)
215
216	c Initialisation pour prendre en compte les vents en 1-D
217	c ------------------------------------------------------
218
219	c vent geostrophique
220	PRINT *,'composante vers l est du vent geostrophique (U) ?'
221	READ(unit,*) gru
222	PRINT *,'composante vers le nord du vent geostrophique (V) ?'
223	READ(unit,*) grv
224
225	c Initialisation des vents au premier pas de temps
226	DO ilayer=1,nlayer
227	u(ilayer)=gru
228	v(ilayer)=grv
229	ENDDO
230
231	c calcul des pressions et altitudes en utilisant les niveaux sigma
232	c ----------------------------------------------------------------
233
234	c Vertical Coordinates (hybrids)
235	c """"""""""""""""""""
236	CALL disvert_noterre
237
238	c Calcul au milieu des couches : Vient de la version Mars
239	c WARNING : le choix de placer le milieu des couches au niveau de
240	c pression intermédiaire est arbitraire et pourrait etre modifié.
241	c C'est fait de la meme facon dans disvert
242
243	DO l = 1, llm
244	aps(l) = 0.5 *( ap(l) +ap(l+1))
245	bps(l) = 0.5 *( bp(l) +bp(l+1))
246	ENDDO
247
248	DO ilevel=1,nlevel
249	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
250	ENDDO
251
252	DO ilayer=1,nlayer
253	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
254	pk(ilayer) =cpp(play(ilayer)/preff)*rcp
255	c write(120,*) ilayer,plev(ilayer),play(ilayer)
256	ENDDO
257	c write(120,*) nlevel,plev(nlevel)
258	c stop
259
260	pks=cpp(psurf/preff)*rcp
261
262	c profil de temperature et altitude au premier appel
263	c --------------------------------------------------
264
265	c modif par rapport a Mars:
266	c on envoie dz/T=-log(play/psurf)*r/g dans profile
267	tmp1(0)=0.0
268	tmp1(1)= -log(play(1)/psurf)*r/g
269	DO ilayer=2,nlayer
270	tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
271	ENDDO
272	call profile(unit,nlayer+1,tmp1,tmp2)
273	CLOSE(unit)
274
275	print*," Pression Altitude Temperature"
276	ilayer=1
277	tsurf(1)=tmp2(0)
278	temp(1)=tmp2(1)
279	zlay(1)=tmp2(1)*tmp1(1)
280	print*," 0",tsurf(1)
281	print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
282	DO ilayer=2,nlayer
283	temp(ilayer)=tmp2(ilayer)
284	zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
285	print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
286	ENDDO
287
288	c temperature du sous-sol
289	c ~~~~~~~~~~~~~~~~~~~~~~~
290	DO isoil=1,nsoil
291	tsoil(isoil)=tsurf(1)
292	ENDDO
293
294	c Initialisation des traceurs
295	c ---------------------------
296
297	DO iq=1,nqtot
298	DO ilayer=1,nlayer
299	q(ilayer,iq) = 0.
300	ENDDO
301	ENDDO
302
303	c Initialisation des parametres d'oro
304	c -----------------------------------
305
306	zmea(1) = 0.
307	zstd(1) = 0.
308	zsig(1) = 0.
309	zgam(1) = 0.
310	zthe(1) = 0.
311	zpic(1) = 0.
312	zval(1) = 0.
313
314	c Ecriture de "startphy.nc"
315	c -------------------------
316	c (Ce fichier sera aussitot relu au premier
317	c appel de "physiq", mais il est necessaire pour passer
318	c les variables purement physiques a "physiq"...
319
320	solsw(1) = 0.
321	sollwdown(1)= 0.
322	dlw(1) = 0.
323	radsol(1) = 0.
324
325	radpas = NINT(1.*day_step/nbapp_rad)
326	soil_model = .true.
327
328	call phyredem("startphy.nc ",
329	. lati,long,
330	. tsurf,tsoil,albedo,
331	. solsw,sollwdown,dlw,radsol,
332	. zmea, zstd, zsig, zgam, zthe, zpic, zval,
333	. temp)
334
335	c=======================================================================
336	c BOUCLE TEMPORELLE DU MODELE 1D
337	c=======================================================================
338	c
339	firstcall=.true.
340	lastcall=.false.
341
342	DO idt=1,ndt
343	IF (idt.eq.ndt) then
344	lastcall=.true.
345	c toujours nulle dans le cas de Venus, pour l'instant...
346	zls = 0.0
347	c write(103,) 'Ls=',zls180./pi
348	c write(103,*) 'Lat=', lati(1)
349	c write(103,*) 'RunEnd - Atmos. Temp. File'
350	c write(103,*) 'RunEnd - Atmos. Temp. File'
351	c write(104,) 'Ls=',zls180./pi
352	c write(104,*) 'Lat=', lati(1)
353	c write(104,*) 'RunEnd - Atmos. Temp. File'
354	ENDIF
355
356	c calcul du geopotentiel
357	c ~~~~~~~~~~~~~~~~~~~~~
358	! ADAPTATION GCM POUR CP(T)
359	DO ilayer=1,nlayer
360	s(ilayer)=(play(ilayer)/psurf)**rcp
361	ENDDO
362	phi(1)=cpptemp(1)(1.E+0-s(1))
363	DO ilayer=2,nlayer
364	phi(ilayer)=phi(ilayer-1)+
365	& cpp(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))0.5
366	& *(s(ilayer-1)-s(ilayer))
367
368	ENDDO
369
370	c appel de la physique
371	c --------------------
372
373	CALL physiq (1,llm,nqtot,
374	, firstcall,lastcall,
375	, day,time,dtphys,
376	, plev,play,pk,phi,phisfi,
377	, presnivs,
378	, u,v,temp,q,
379	, w,
380	C - sorties
381	s du,dv,dtemp,dq,dpsurf)
382
383	c print*,"DT APRES PHYSIQ=",day,time
384	c print*,dtemp
385	c print*,temp
386	c print*," "
387	c stop
388
389	c evolution du vent : modele 1D
390	c -----------------------------
391
392	c la physique calcule les derivees temporelles de u et v.
393	c Pas de coriolis
394	DO ilayer=1,nlayer
395	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
396	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
397	ENDDO
398	c
399	c Calcul du temps au pas de temps suivant
400	c ---------------------------------------
401	firstcall=.false.
402	time=time+dtphys/daysec
403	IF (time.gt.1.E+0) then
404	time=time-1.E+0
405	day=day+1
406	ENDIF
407
408	c calcul des vitesses et temperature au pas de temps suivant
409	c ----------------------------------------------------------
410
411	DO ilayer=1,nlayer
412	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
413	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
414	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
415	ENDDO
416
417	c calcul des pressions au pas de temps suivant
418	c ----------------------------------------------------------
419
420	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
421	DO ilevel=1,nlevel
422	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
423	ENDDO
424	DO ilayer=1,nlayer
425	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
426	ENDDO
427
428	ENDDO ! fin de la boucle temporelle
429
430	c ========================================================
431	c GESTION DES SORTIE
432	c ========================================================
433
434	print*,"Temperature finale:"
435	print*,temp
436
437	c stabilite
438	DO ilayer=1,nlayer
439	zlay(ilayer) = phi(ilayer)/g/1000. !en km
440	ENDDO
441	DO ilayer=2,nlayer
442	tmp1(ilayer) =
443	. (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1))
444	. + 1000.*g/cpp
445	ENDDO
446
447	OPEN(11,file='profile.new')
448	write (11,*) tsurf
449	DO ilayer=1,nlayer
450	write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
451	ENDDO
452
453	c ========================================================
454	END
455
456	c***********************************************************************
457	c***********************************************************************
458
459	#include "../dyn3d/disvert_noterre.F"
460	#include "../dyn3d/abort_gcm.F"
461	#include "../dyn3d/dump2d.F"
462	#include "../dyn3d/cpdet.F"
463

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

Download in other formats:

Original Format