Context Navigation

rcm1d.F @ 3094

Last change on this file since 3094 was 2851, checked in by streel, 2 years ago

LMDZ.VENUS - chemical reaction rates update with JPL 2019

Update of the main 1D subroutine with nitrogen species + generalisation of VMR output every 1/5 venusian day
Added the H2 photolysis
Added the traceur.def for all neutral species only
Added the VMR def files for 1D model

File size: 16.4 KB

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

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

Download in other formats:

Original Format