Context Navigation

rcm1d.F @ 1443

Last change on this file since 1443 was 1443, checked in by emillour, 10 years ago
Titan and Venus GCMs: Follow-up to the changes in dynamics/physics interface: ener.h, logic.h, serre.h and temps.h are now modules. EM
File size: 13.6 KB

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

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

Download in other formats:

Original Format