Context Navigation

rcm1d.F @ 3556

Last change on this file since 3556 was 3540, checked in by slebonnois, 2 weeks ago
SL: bug for 1D (was using constant Cp since cpofT keyword not defined)
File size: 16.0 KB

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

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

Download in other formats:

Original Format