Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

testphys1d.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 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.3 KB

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

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

Download in other formats:

Original Format