Context Navigation

testphys1d_triton.F

Last change on this file was 3175, checked in by emillour, 10 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 26.1 KB

Line
1	PROGRAM testphys1d_triton
2
3	! to use 'getin'
4	USE ioipsl_getincom
5	use planet_h
6	use comgeomfi_h
7	use comsoil_h
8	IMPLICIT NONE
9
10	!==================================================================
11	!
12	! Purpose
13	! -------
14	! Run the physics package of the universal model in a 1D column.
15	!
16	! It can be compiled with a command like (e.g. for 25 layers):
17	! "makegcm -p pluto -d 25 testphys1d"
18	! It requires the files "callphys.def", "z2sig.def",
19	! "traceur.def", and "run1d.def" with a line "INCLUDEDEF=callphys.def"
20	!
21	! Authors
22	! -------
23	! Frederic Hourdin
24	! R. Fournier
25	! F. Forget
26	! F. Montmessin (water ice added)
27
28	!==================================================================
29
30	#include "dimensions.h"
31	#include "dimphys.h"
32	#include "surfdat.h"
33	!#include "comsoil.h"
34	#include "comdiurn.h"
35	#include "callkeys.h"
36	#include "comcstfi.h"
37	#include "comsaison.h"
38	#include "control.h"
39	#include "comvert.h"
40	#include "netcdf.inc"
41	#include "comg1d.h"
42	#include "fisice.h"
43	#include "logic.h"
44	#include "advtrac.h"
45
46	c --------------------------------------------------------------
47	c Declarations
48	c --------------------------------------------------------------
49	c
50	INTEGER nlayer,nlevel,nsoil,ndt
51	INTEGER ilayer,ilevel,isoil,idt,iq
52	LOGICAl firstcall,lastcall
53	c
54	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
55	INTEGER lecttsoil ! lecture of tsoil from proftsoil
56	INTEGER lecthaze ! lecture of haze from profhaze
57	REAL day ! date durant le run
58	REAL AU ! astronomical unit AU=149 million km
59	REAL time ! time (0<time<1 ; time=0.5 a midi)
60	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
61	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
62	REAL psurf,tsurf
63	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
64	REAL gru,grv ! prescribed "geostrophic" background wind
65	REAL temp(nlayermx) ! temperature at the middle of the layers
66	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
67	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
68	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
69	integer :: i_n2=0 ! tracer index of n2 ice
70	integer :: i_ch4_ice=0 ! tracer index of Ch4 ice
71	integer :: i_ch4_gas=0 ! tracer index of Ch4 gas
72	integer :: i_co_ice=0 ! tracer index of Co ice
73	integer :: i_co_gas=0 ! tracer index of Co gas
74	integer :: i_prec_haze=0 ! tracer index of haze precursor
75	integer :: i_haze=0 ! tracer index of haze
76	integer :: i_haze10=0 ! tracer index of haze
77	integer :: i_haze30=0 ! tracer index of haze
78	integer :: i_haze50=0 ! tracer index of haze
79	integer :: i_haze100=0 ! tracer index of haze
80	REAL emis ! surface layer
81	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
82	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
83	REAL Lsperi,excentric
84
85	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
86	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
87	REAL dpsurf
88	REAL dq(nlayermx,nqmx)
89
90	c Various intermediate variables
91	REAL zls
92	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
93	REAL pks, ptif, w(nlayermx)
94	INTEGER ierr
95	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
96	Logical tracerdyn
97	integer :: nq=1 ! number of tracers
98
99	character (len=7) :: str7
100
101	logical saveprofile
102
103	! added by RW for zlay computation
104	real Hscale, rho, dz
105
106	real ch4ref ! % ch4
107	real coref ! % co
108	real hazeref ! haze kg/kg
109	real prechazeref ! prec haze kg/kg
110
111
112	c=======================================================================
113
114	c=======================================================================
115	c INITIALISATION
116	c=======================================================================
117
118	saveprofile=.true.
119
120	c ------------------------------------------------------
121	c Default values for constants (corresponding to Triton)
122	c ------------------------------------------------------
123
124	pi=2.E+0*asin(1.E+0)
125
126	c Constante de Triton
127	c ----------------------------
128	rad=1353400. ! rayon de Triton (m)
129	daysec=507773 ! duree du sol (s)
130	omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1)
131	g=0.78 ! gravite (m.s-2)
132	mugaz=28. ! Masse molaire de l'atm (g.mol-1) ~28
133	rcp=.2853 ! = r/cp ~0.2853
134	r= 8.314511E+0 *1000.E+0/mugaz
135	cpp= r/rcp
136	year_day = 10247 ! duree de l'annee (sols)
137	periheli = 4444.45 ! dist.min. soleil-Triton (Mkm)
138	aphelie = 4545.67 ! dist.max. soleil-Triton (Mkm)
139	peri_day = 0 ! date du perihelie (sols depuis printemps)
140	obliquit=28.32 ! Obliquite de la planete (deg)
141	excentric= 0. ! Excentricite
142
143	c Parametres Couche limite et Turbulence
144	c --------------------------------------
145	z0 = 1.e-2 ! surface roughness (m) ~0.01
146
147	c ----------------------------------------------------
148	hybrid=.true.
149
150	! Constants for stokes.F90
151	molrad=1.93e-10 ! N2 TB16!!!
152	visc=6.67e-6 ! N2 TB16!!!
153
154	c ------------------------------------------------------
155	c Load parameters from "run.def"
156	c ------------------------------------------------------
157	open(90,file='run.def',status='old',form='formatted',
158	& iostat=ierr)
159	if (ierr.ne.0) then
160	write(,) 'Cannot find required file "run.def"'
161	write(,) ' (which should contain some input parameters'
162	write(,) ' along with the following line:'
163	write(,) ' INCLUDEDEF=callphys.def'
164	write(,) ' )'
165	write(,) ' ... might as well stop here ...'
166	stop
167	else
168	close(90)
169	endif
170
171	! check if we are going to run with or without tracers
172	write(,) "Run with or without tracer transport ?"
173	tracer=.false. ! default value
174	call getin("tracer",tracer)
175	write(,) " tracer = ",tracer
176
177	! while we're at it, check if there is a 'traceur.def' file
178	! and preocess it, if necessary. Otherwise initialize tracer names
179	if (tracer) then
180	! load tracer names from file 'traceur.def'
181	open(90,file='traceur.def',status='old',form='formatted',
182	& iostat=ierr)
183	if (ierr.eq.0) then
184	write(,) "testphys1d: Reading file traceur.def"
185	! read number of tracers:
186	read(90,*,iostat=ierr) nq
187	if (ierr.ne.0) then
188	write(,) "testphys1d: error reading number of tracers"
189	write(,) " (first line of traceur.def) "
190	stop
191	else
192	! check that the number of tracers is indeed nqmx
193	if (nq.ne.nqmx) then
194	write(,) "testphys1d: error, wrong number of tracers:"
195	write(,) "nq=",nq," whereas nqmx=",nqmx
196	stop
197	endif
198	endif
199
200	! initialize advection schemes to Van-Leer for all tracers
201	do iq=1,nq
202	iadv(iq)=3 ! Van-Leer
203	enddo
204
205	do iq=1,nq
206	! minimal version, just read in the tracer names, 1 per line
207	read(90,*,iostat=ierr) tnom(iq)
208	write(,) 'Ini 1d reading traceur.def: ', tnom(iq),iq
209	if (ierr.ne.0) then
210	write(,) 'testphys1d: error reading tracer names...'
211	stop
212	endif
213	enddo !of do iq=1,nq
214
215	! TB check for n2_ice / ch4 tracers:
216	write(,) 'Tracers: tnom=', tnom(:)
217	do iq=1,nq
218	if (tnom(iq)=="n2") then
219	i_n2=iq
220	elseif (tnom(iq)=="ch4_ice") then
221	i_ch4_ice=iq
222	elseif (tnom(iq)=="ch4_gas") then
223	i_ch4_gas=iq
224	elseif (tnom(iq)=="prec_haze") then
225	i_prec_haze=iq
226	elseif (tnom(iq)=="haze") then
227	i_haze=iq
228	elseif (tnom(iq)=="haze10") then
229	i_haze10=iq
230	elseif (tnom(iq)=="haze30") then
231	i_haze30=iq
232	elseif (tnom(iq)=="haze50") then
233	i_haze50=iq
234	elseif (tnom(iq)=="haze100") then
235	i_haze100=iq
236	elseif (tnom(iq)=="co_ice") then
237	i_co_ice=iq
238	elseif (tnom(iq)=="co_gas") then
239	i_co_gas=iq
240	endif
241	enddo
242
243	else ! ierr=0
244	write(,) 'Cannot find required file "traceur.def"'
245	write(,) ' If you want to run with tracers, I need it'
246	write(,) ' ... might as well stop here ...'
247	stop
248	endif
249	close(90)
250
251	else
252	! we still need to set (dummy) tracer names for physdem1
253	nq=nqmx
254	do iq=1,nq
255	write(str7,'(a1,i2.2)')'q',iq
256	tnom(iq)=str7
257	enddo
258	! actually, we'll need at least one "n2_ice" tracer
259	! (for surface N2 ice)
260	write(,) "No tracer ! tracer=false"
261	tnom(1)="n2"
262	i_n2=1
263	endif ! of if (tracer)
264
265	c Date et heure locale du debut du run
266	c ------------------------------------
267	c Date (en sols depuis le solstice de printemps) du debut du run
268	day0 = 0 ! default value for day0
269	write(,) 'Initial date (in sols ; =0 at Ls=0)?'
270	call getin("day0",day0)
271	day=float(day0)
272	write(,) " day0 = ",day0
273	c Heure de demarrage
274	time=0 ! default value for time
275	write(,)'Initial local time (in hours, between 0 and 24)?'
276	call getin("time",time)
277	write(,)" time = ",time
278	time=time/24.E+0 ! convert time (hours) to fraction of sol
279
280	c Discretisation (Definition de la grille et des pas de temps)
281	c --------------
282	c
283	nlayer=nlayermx
284	nlevel=nlayer+1
285	nsoil=nsoilmx
286	PRINT *,'Dims nlevel=',nlevel,' nsoil=',nsoil
287
288	PRINT *,'Length of day (s) ?'
289	call getin("daysec",daysec)
290	write(,) " daysec = ",daysec
291	omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1)
292
293	PRINT *,'Length of year (sol) ?'
294	call getin("year_day",year_day)
295	write(,) " year_day= ",year_day
296
297	day_step=48 ! default value for day_step
298	PRINT *,'Number of time steps per sol ?'
299	call getin("day_step",day_step)
300	write(,) " day_step = ",day_step
301
302	ndt=10 ! default value for ndt
303	PRINT *,'Number of sols to run ?'
304	call getin("ndt",ndt)
305	write(,) " ndt = ",ndt
306
307	ndt=ndt*day_step
308	dtphys=daysec/day_step
309
310	c Pression de surface sur la planete
311	c ------------------------------------
312	psurf=1. ! default value for psurf
313	PRINT *,'Surface pressure (Pa) ?'
314	call getin("psurf",psurf)
315	write(,) " psurf = ",psurf
316	c Pression de reference
317	preff=1. ! these values are not needed in 1D
318	pa=0.25*preff
319	c Gravity of planet
320	c -----------------
321	g=0.78 ! default value for g
322	PRINT *,'Gravity ?'
323	call getin("g",g)
324	write(,) " g = ",g
325
326	c Molar mass of atmosphere
327	c ------------------------
328	PRINT *,'Molar mass of atmosphere ?'
329	call getin("mugaz",mugaz)
330	write(,) " mugaz = ",mugaz
331
332	c Specific heat capacity of atmosphere
333	c ------------------------------------
334	PRINT *,'Specific heat capacity of atmosphere ?'
335	call getin("cpp",cpp)
336	write(,) " cpp = ",cpp
337
338	r= 8.314511E+0 *1000.E+0/mugaz
339	rcp=r/cpp
340
341	c latitude/longitude
342	c ------------------
343	lati(1)=0 ! default value for lati(1)
344	PRINT *,'latitude (in degrees) ?'
345	call getin("latitude",lati(1))
346	write(,) " latitude = ",lati(1)
347	lati(1)=lati(1)*pi/180.E+0
348	long(1)=0.E+0
349	long(1)=long(1)*pi/180.E+0
350
351	c periastron/apoastron
352	c --------------------
353	PRINT *,'periastron (in AU) ?'
354	call getin("periheli",periheli)
355	write(,) "periastron = ",periheli
356	AU=149.597870700 ! km
357	periheli=periheli*AU ! AU to Mkm
358	PRINT *,'apoastron (in AU) ?'
359	call getin("aphelie",aphelie)
360	write(,) "apoastron = ",aphelie
361	aphelie=aphelie*AU ! AU to Mkm
362
363	excentric=(1-periheliperiheli/(aphelieaphelie) )**0.5
364
365	c obliquity
366	c ---------
367	PRINT *,'obliquity (in deg) ?'
368	call getin("obliquit",obliquit)
369	write(,) "obliquity = ",obliquit
370
371	c Initialisation albedo /emis / inertie du sol
372	c --------------------------------------
373	albedodat(1)=0.2 ! default value for albedodat
374	PRINT *,'Albedo of bare ground ?'
375	call getin("albedo",albedodat(1))
376	write(,) " albedo = ",albedodat(1)
377
378	emis=0.8 ! default value for emissivity
379	PRINT *,'Emissivity of bare ground ?'
380	call getin("emis",emis)
381	write(,) " emis = ",emis
382
383	inertiedat(1,1)=400 ! default value for inertiedat
384	PRINT *,'Soil thermal inertia (SI) ?'
385	call getin("inertia",inertiedat(1,1))
386	write(,) " inertia = ",inertiedat(1,1)
387	c
388	c pour le schema d'ondes de gravite
389	c ---------------------------------
390	zmea(1)=0.E+0
391	zstd(1)=0.E+0
392	zsig(1)=0.E+0
393	zgam(1)=0.E+0
394	zthe(1)=0.E+0
395
396	c Initialisation des traceurs
397	c ---------------------------
398	DO iq=1,nqmx
399	if (iq.eq.i_n2) then
400	DO ilayer=1,nlayer
401	q(ilayer,iq) = 1
402	ENDDO
403	write(,) 'ini 1D traceur ',iq,' : q_n2 = ', q(:,iq)
404	else if (iq.eq.i_ch4_gas) then
405	ch4ref=0.5 ! default value for vmr ch4
406	PRINT *,'vmr CH4 (%) ?'
407	call getin("ch4ref",ch4ref)
408	write(,) " CH4 ref (%) = ",ch4ref
409	DO ilayer=1,nlayer
410	q(ilayer,iq) = 0.01ch4ref16./28.
411	ENDDO
412	write(,) 'ini 1D traceur ',iq,' : q_ch4_gas = ', q(:,iq)
413	else if (iq.eq.i_co_gas) then
414	coref=0.05 ! default value for vmr ch4
415	PRINT *,'vmr CO (%) ?'
416	call getin("coref",coref)
417	write(,) " CO ref (%) = ",coref
418	DO ilayer=1,nlayer
419	q(ilayer,iq) = 0.01coref28./28.
420	ENDDO
421	write(,) 'ini 1D traceur ',iq,' : q_co_gas = ', q(:,iq)
422	else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.
423	& (iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then
424	hazeref=0. ! default value for haze mmr
425	PRINT *,'qhaze (kg/kg) ?'
426	call getin("hazeref",hazeref)
427	write(,) " haze ref (kg/kg) = ",hazeref
428	DO ilayer=1,nlayer
429	q(ilayer,iq) = hazeref
430	ENDDO
431	write(,) 'ini 1D traceur ',iq,' : q_haze = ', q(:,iq)
432	else if (iq.eq.i_prec_haze) then
433	prechazeref=0. ! default value for vmr ch4
434	PRINT *,'qprechaze (kg/kg) ?'
435	call getin("prechazeref",prechazeref)
436	write(,) " prec haze ref (kg/kg) = ",prechazeref
437	DO ilayer=1,nlayer
438	q(ilayer,iq) = prechazeref
439	ENDDO
440	write(,) 'ini 1D traceur ',iq,' : q_prec_haze = ', q(:,iq)
441
442	else
443	DO ilayer=1,nlayer
444	q(ilayer,iq) = 0.
445	ENDDO
446	endif
447	ENDDO
448
449	c TB: clean surface
450	DO iq=1,nqmx
451	qsurf(iq) = 0.
452	ENDDO
453
454	c Initialisation speciales "physiq"
455	c ---------------------------------
456	c la surface de chaque maille est inutile en 1D --->
457	area(1)=1.E+0
458
459	c le geopotentiel au sol est inutile en 1D car tout est controle
460	c par la pression de surface --->
461	phisfi(1)=0.E+0
462
463	c "inifis" reproduit un certain nombre d'initialisations deja faites
464	c + lecture des clefs de callphys.def
465	c + calcul de la frequence d'appel au rayonnement
466	c + calcul des sinus et cosinus des longitude latitude
467
468	! Possible issue with dtphys in input and include!!!
469	CALL inifis(1,llm,day0,daysec,dtphys,
470	. lati,long,area,rad,g,r,cpp)
471	c Initialisation pour prendre en compte les vents en 1-D
472	c ------------------------------------------------------
473	ptif=2.E+0omegsinlat(1)
474
475	c vent geostrophique
476	gru=10. ! default value for gru
477	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
478	call getin("u",gru)
479	write(,) " u = ",gru
480	grv=0. !default value for grv
481	PRINT *,'meridional northward component of the geostrophic',
482	&' wind (m/s) ?'
483	call getin("v",grv)
484	write(,) " v = ",grv
485
486	c Initialisation des vents au premier pas de temps
487	DO ilayer=1,nlayer
488	u(ilayer)=gru
489	v(ilayer)=grv
490	ENDDO
491
492	c energie cinetique turbulente
493	DO ilevel=1,nlevel
494	q2(ilevel)=0.E+0
495	ENDDO
496
497	c Surface
498	c -------------------------------------------
499	if(i_n2.gt.0)then
500	qsurf(i_n2)=0 ! default value for N2ice
501	print*,'Initial N2 ice on the surface (kg.m-2)'
502	call getin("n2ice",qsurf(i_n2))
503	write(,) " n2ice = ",qsurf(i_n2)
504	endif
505
506	c calcul des pressions et altitudes en utilisant les niveaux sigma
507	c ----------------------------------------------------------------
508
509	c Vertical Coordinates
510	c """"""""""""""""""""
511	hybrid=.false.
512	PRINT *,'Hybrid coordinates ?'
513	call getin("hybrid",hybrid)
514	write(,) " hybrid = ", hybrid
515
516	CALL disvert
517
518	! we want only the scale height from z2sig, in order to compute zlay
519	open(91,file="z2sig.def",form='formatted')
520	read(91,*) Hscale
521	close(91)
522
523	DO ilevel=1,nlevel
524	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
525	ENDDO
526
527	DO ilayer=1,nlayer
528	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
529	ENDDO
530
531	DO ilayer=1,nlayer
532	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
533	! & /g
534	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
535	ENDDO
536
537	c Orbital parameters
538	c ------------------
539	PRINT *,'Solar longitude of perihelion Lsperi '
540	call getin("Lsperi",Lsperi)
541	write(,) " Lsperi = ", Lsperi
542
543	Lsperi = Lsperi * pi/180. ! Ls of perihelion
544
545	c Calcul du sol correspondant a Lsperi
546	call call_dayperi(Lsperi,excentric,peri_day,year_day)
547	PRINT *,'date of perihelion (sol)',peri_day
548
549	c Profil de temperature au premier appel
550	c --------------------------------------
551	pks=psurf**rcp
552
553	c Altitude en km dans profile: on divise zlay par 1000
554	tmp1(0)=0.E+0
555	DO ilayer=1,nlayer
556	tmp1(ilayer)=zlay(ilayer)/1000.E+0
557	ENDDO
558	call profile(nlayer+1,tmp1,tmp2)
559
560	tsurf=tmp2(0)
561	DO ilayer=1,nlayer
562	temp(ilayer)=tmp2(ilayer)
563	ENDDO
564
565	! Initialize soil properties and temperature
566	! ------------------------------------------
567	volcapa=1.e6 ! volumetric heat capacity
568
569	lecttsoil=0 ! default value for lecttsoil
570	call getin("lecttsoil",lecttsoil)
571	if (lecttsoil==1) then
572	OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
573	DO isoil=1,nsoil
574	READ (14,*) tsoil(isoil)
575	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
576	ENDDO
577	GOTO 201
578	101 STOP'fichier proftsoil inexistant'
579	201 CONTINUE
580	CLOSE(14)
581
582	else
583	DO isoil=1,nsoil
584	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
585	tsoil(isoil)=tsurf ! soil temperature
586	ENDDO
587	endif
588
589	! Initialize depths
590	! -----------------
591	do isoil=0,nsoil-1
592	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
593	enddo
594	do isoil=1,nsoil
595	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
596	enddo
597
598	c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
599	c ----------------------------------------------------------------
600	c (effectuee avec "writeg1d" a partir de "physiq.F"
601	c ou tout autre subroutine physique, y compris celle ci).
602
603	g1d_nlayer=nlayer
604	g1d_nomfich='g1d.dat'
605	g1d_unitfich=40
606	g1d_nomctl='g1d.ctl'
607	g1d_unitctl=41
608	g1d_premier=.true.
609	g2d_premier=.true.
610
611	! Initialize haze profile
612	! ------------------------------------------
613	if (haze) then
614
615	lecthaze=0 ! default value for lecthaze
616	call getin("lecthaze",lecthaze)
617	if (lecthaze==1) then
618	OPEN(15,file='profhaze',status='old',form='formatted',err=301)
619	DO iq=1,nq
620	if (iq.eq.i_haze) then
621	DO ilayer=1,nlayer
622	READ (15,*) q(ilayer,iq)
623	ENDDO
624	endif
625	ENDDO
626	GOTO 401
627	301 STOP'fichier profhaze inexistant'
628	401 CONTINUE
629	CLOSE(15)
630	endif
631	endif
632
633	c Ecriture de "startfi"
634	c --------------------
635	c (Ce fichier sera aussitot relu au premier
636	c appel de "physiq", mais il est necessaire pour passer
637	c les variables purement physiques a "physiq"...
638
639	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
640	. dtphys,float(day0),
641	. time,tsurf,tsoil,emis,q2,qsurf,
642	. area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
643
644	c=======================================================================
645	c=======================================================================
646	c=======================================================================
647	c BOUCLE TEMPORELLE DU MODELE 1D
648	c=======================================================================
649	c=======================================================================
650	c=======================================================================
651
652	firstcall=.true.
653	lastcall=.false.
654
655	DO idt=1,ndt
656	IF (idt.eq.ndt-day_step-1) then !test
657	lastcall=.true.
658	call solarlong(day*1.0,zls)
659	write(103,) 'Ls=',zls180./pi
660	write(103,) 'Lat=', lati(1)180./pi
661	write(103,*) 'RunEnd - Atmos. Temp. File'
662	write(103,*) 'RunEnd - Atmos. Temp. File'
663	write(104,) 'Ls=',zls180./pi
664	write(104,*) 'Lat=', lati(1)
665	write(104,*) 'RunEnd - Atmos. Temp. File'
666	ENDIF
667
668	c calcul du geopotentiel
669	c ~~~~~~~~~~~~~~~~~~~~~
670	DO ilayer=1,nlayer
671	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
672	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
673	ENDDO
674
675	phi(1)=pksh(1)(1.E+0-s(1))
676	DO ilayer=2,nlayer
677	phi(ilayer)=phi(ilayer-1)+
678	& pks(h(ilayer-1)+h(ilayer)).5E+0
679	& *(s(ilayer-1)-s(ilayer))
680	ENDDO
681
682	c appel de la physique
683	c --------------------
684	CALL physiq (1,llm,nqmx,
685	, firstcall,lastcall,
686	, day,time,dtphys,
687	, plev,play,phi,
688	, u, v,temp, q,
689	, w,
690	C - sorties
691	s du, dv, dtemp, dq,dpsurf,tracerdyn)
692
693
694	c evolution du vent : modele 1D
695	c -----------------------------
696
697	c la physique calcule les derivees temporelles de u et v.
698	c on y rajoute betement un effet Coriolis.
699	c
700	c DO ilayer=1,nlayer
701	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
702	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
703	c ENDDO
704
705	c Pour certain test : pas de coriolis a l'equateur
706	c if(lati(1).eq.0.) then
707	DO ilayer=1,nlayer
708	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
709	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
710	ENDDO
711	c end if
712	c
713	c Calcul du temps au pas de temps suivant
714	c ---------------------------------------
715	firstcall=.false.
716	time=time+dtphys/daysec
717	IF (time.gt.1.E+0) then
718	time=time-1.E+0
719	day=day+1
720	ENDIF
721
722	c calcul des vitesses et temperature au pas de temps suivant
723	c ----------------------------------------------------------
724
725	DO ilayer=1,nlayer
726	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
727	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
728	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
729	ENDDO
730
731	c calcul des pressions au pas de temps suivant
732	c ----------------------------------------------------------
733
734	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
735	DO ilevel=1,nlevel
736	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
737	ENDDO
738	DO ilayer=1,nlayer
739	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
740	ENDDO
741
742	c calcul traceur au pas de temps suivant
743	c --------------------------------------
744	DO iq = 1, nqmx
745	DO ilayer=1,nlayer
746	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
747	ENDDO
748	END DO
749
750	ENDDO ! fin de la boucle temporelle
751
752	c ========================================================
753	c GESTION DES SORTIE
754	c ========================================================
755
756	c fermeture pour conclure les sorties format grads dans "g1d.dat"
757	c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
758	c ou tout autre subroutine physique, y compris celle ci).
759
760	! save atm temperature profile
761	if(saveprofile)then
762	OPEN(12,file='profile.out',form='formatted')
763	write(12,*) tsurf
764	DO ilayer=1,nlayermx
765	write(12,*) temp(ilayer)
766	ENDDO
767	CLOSE(12)
768	endif
769	! save haze profile
770	if (haze.and.lecthaze.eq.1) then
771	OPEN(16,file='profhaze.out',form='formatted')
772	DO iq=1,nq
773	if (iq.eq.i_haze) then
774	DO ilayer=1,nlayer
775	write(16,*) q(ilayer,iq)
776	ENDDO
777	endif
778	ENDDO
779	CLOSE(16)
780	endif
781
782	! here we compute altitude CORRECTLY!!!
783
784	! print*,'zlay=',zlay/1000.
785	! rho = play(1)/(R*tsurf)
786	! zlay(1)=(plev(1)-plev(2))/(g*rho)
787	! do ilayer=2,nlayer
788	! rho = play(ilayer)/(R*temp(ilayer))
789	! rho = play(ilayer)/(R*230.0)
790	! dz = (plev(ilayer-1)-plev(ilayer))/(g*rho)
791	! zlay(ilayer)=zlay(ilayer-1)+dz
792	! enddo
793	! print*,'zlay=',zlay/1000.
794	! stop
795
796	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
797	do ilayer=1,nlayer
798	zlay(ilayer)=-18.0e3*log(play(ilayer)/psurf)
799	enddo
800	CALL endg1d(1,nlayer,zlay/1000.,ndt)
801
802	c ========================================================
803	END
804
805	c***********************************************************************
806	c***********************************************************************
807	c Subroutines Bidons utilise seulement en 3D, mais
808	c necessaire a la compilation de testphys1d en 1-D
809
810	subroutine gr_fi_dyn
811	RETURN
812	END
813
814	c***********************************************************************
815	c***********************************************************************
816
817	#include "../dyn3d/disvert.F"
818	#include "call_dayperi.F"

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

Download in other formats:

Original Format