Context Navigation

rcm1d.F @ 1395

Last change on this file since 1395 was 1395, checked in by emillour, 10 years ago
All GCMS: Some cleanup and tidying on the dynamics/physics interface. Essentially affects the "iniphysiq" routine in all physics packages. EM
File size: 13.7 KB

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

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

Download in other formats:

Original Format