Context Navigation

testphys1d.F @ 706

Last change on this file since 706 was 627, checked in by emillour, 13 years ago

Mars GCM:

Some cleanup on messages and comments in code about the reference pressure

for dust opacity which is now 610 Pa.

File size: 24.9 KB

Line
1
2	PROGRAM testphys1d
3	! to use 'getin'
4	USE ioipsl_getincom
5	IMPLICIT NONE
6
7	c=======================================================================
8	c subject:
9	c --------
10	c PROGRAM useful to run physical part of the martian GCM in a 1D column
11	c
12	c Can be compiled with a command like (e.g. for 25 layers)
13	c "makegcm -p mars -d 25 testphys1d"
14	c It requires the files "testphys1d.def" "callphys.def"
15	c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
16	c and a file describing the sigma layers (e.g. "z2sig.def")
17	c
18	c author: Frederic Hourdin, R.Fournier,F.Forget
19	c -------
20	c
21	c update: 12/06/2003 including chemistry (S. Lebonnois)
22	c and water ice (F. Montmessin)
23	c
24	c=======================================================================
25
26	#include "dimensions.h"
27	#include "dimphys.h"
28	#include "dimradmars.h"
29	#include "comgeomfi.h"
30	#include "surfdat.h"
31	#include "comsoil.h"
32	#include "comdiurn.h"
33	#include "callkeys.h"
34	#include "comcstfi.h"
35	#include "planete.h"
36	#include "comsaison.h"
37	#include "yomaer.h"
38	#include "control.h"
39	#include "comvert.h"
40	#include "netcdf.inc"
41	#include "comg1d.h"
42	#include "logic.h"
43	#include "advtrac.h"
44
45	c --------------------------------------------------------------
46	c Declarations
47	c --------------------------------------------------------------
48	c
49	INTEGER unitstart ! unite d'ecriture de "startfi"
50	INTEGER nlayer,nlevel,nsoil,ndt
51	INTEGER ilayer,ilevel,isoil,idt,iq
52	LOGICAl firstcall,lastcall
53	c
54	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
55	c
56	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
57	REAL day ! date durant le run
58	REAL time ! time (0<time<1 ; time=0.5 a midi)
59	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
60	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
61	REAL psurf,tsurf
62	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
63	REAL gru,grv ! prescribed "geostrophic" background wind
64	REAL temp(nlayermx) ! temperature at the middle of the layers
65	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
66	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
67	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
68	REAL co2ice ! co2ice layer (kg.m-2)
69	REAL emis ! surface layer
70	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
71	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
72
73	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
74	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
75	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
76	REAL dpsurf
77	REAL dq(nlayermx,nqmx)
78	REAL dqdyn(nlayermx,nqmx)
79
80	c Various intermediate variables
81	INTEGER thermo
82	REAL zls
83	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
84	REAL pks, ptif, w(nlayermx)
85	REAL qtotinit, mqtot(nqmx),qtot
86	INTEGER ierr, aslun
87	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
88	Logical tracerdyn
89	integer :: nq=1 ! number of tracers
90
91	character*2 str2
92	character (len=7) :: str7
93	character(len=44) :: txt
94
95	c=======================================================================
96
97	c=======================================================================
98	c INITIALISATION
99	c=======================================================================
100
101	c ------------------------------------------------------
102	c Constantes prescrites ICI
103	c ------------------------------------------------------
104
105	pi=2.E+0*asin(1.E+0)
106
107	c Constante de la Planete Mars
108	c ----------------------------
109	rad=3397200. ! rayon de mars (m) ~3397200 m
110	daysec=88775. ! duree du sol (s) ~88775 s
111	omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1)
112	g=3.72 ! gravite (m.s-2) ~3.72
113	mugaz=43.49 ! Masse molaire de l'atm (g.mol-1) ~43.49
114	rcp=.256793 ! = r/cp ~0.256793
115	r= 8.314511E+0 *1000.E+0/mugaz
116	cpp= r/rcp
117	year_day = 669 ! duree de l'annee (sols) ~668.6
118	periheli = 206.66 ! dist.min. soleil-mars (Mkm) ~206.66
119	aphelie = 249.22 ! dist.max. soleil-mars (Mkm) ~249.22
120	peri_day = 485. ! date du perihelie (sols depuis printemps)
121	obliquit = 25.2 ! Obliquite de la planete (deg) ~25.2
122
123	c Parametres Couche limite et Turbulence
124	c --------------------------------------
125	z0_default = 1.e-2 ! surface roughness (m) ~0.01
126	emin_turb = 1.e-6 ! energie minimale ~1.e-8
127	lmixmin = 30 ! longueur de melange ~100
128
129	c propriete optiques des calottes et emissivite du sol
130	c ----------------------------------------------------
131	emissiv= 0.95 ! Emissivite du sol martien ~.95
132	emisice(1)=0.95 ! Emissivite calotte nord
133	emisice(2)=0.95 ! Emissivite calotte sud
134	albedice(1)=0.5 ! Albedo calotte nord
135	albedice(2)=0.5 ! Albedo calotte sud
136	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
137	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
138	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
139	dtemisice(2) = 2. ! time scale for snow metamorphism (south
140	hybrid=.false.
141
142	c ------------------------------------------------------
143	c Lecture des parametres dans "run.def"
144	c ------------------------------------------------------
145
146
147	! check if 'run.def' file is around (otherwise reading parameters
148	! from callphys.def via getin() routine won't work.
149	open(99,file='run.def',status='old',form='formatted',
150	& iostat=ierr)
151	if (ierr.ne.0) then
152	write(,) 'Cannot find required file "run.def"'
153	write(,) ' (which should contain some input parameters'
154	write(,) ' along with the following line:'
155	write(,) ' INCLUDEDEF=callphys.def'
156	write(,) ' )'
157	write(,) ' ... might as well stop here ...'
158	stop
159	else
160	close(99)
161	endif
162
163	! check if we are going to run with or without tracers
164	write(,) "Run with or without tracer transport ?"
165	tracer=.false. ! default value
166	call getin("tracer",tracer)
167	write(,) " tracer = ",tracer
168
169	! while we're at it, check if there is a 'traceur.def' file
170	! and preocess it, if necessary. Otherwise initialize tracer names
171	if (tracer) then
172	! load tracer names from file 'traceur.def'
173	open(90,file='traceur.def',status='old',form='formatted',
174	& iostat=ierr)
175	if (ierr.ne.0) then
176	write(,) 'Cannot find required file "traceur.def"'
177	write(,) ' If you want to run with tracers, I need it'
178	write(,) ' ... might as well stop here ...'
179	stop
180	else
181	write(,) "testphys1d: Reading file traceur.def"
182	! read number of tracers:
183	read(90,*,iostat=ierr) nq
184	if (ierr.ne.0) then
185	write(,) "testphys1d: error reading number of tracers"
186	write(,) " (first line of traceur.def) "
187	stop
188	else
189	! check that the number of tracers is indeed nqmx
190	if (nq.ne.nqmx) then
191	write(,) "testphys1d: error, wrong number of tracers:"
192	write(,) "nq=",nq," whereas nqmx=",nqmx
193	stop
194	endif
195	endif
196	endif
197	! read tracer names from file traceur.def
198	do iq=1,nqmx
199	read(90,*,iostat=ierr) tnom(iq)
200	if (ierr.ne.0) then
201	write(,) 'testphys1d: error reading tracer names...'
202	stop
203	endif
204	enddo
205	close(90)
206
207	! initialize tracers here:
208	write(,) "testphys1d: initializing tracers"
209	q(:,:)=0 ! default, set everything to zero
210	qsurf(:)=0
211	! "smarter" initialization of some tracers
212	! (get values from "profile_*" files, if these are available)
213	do iq=1,nqmx
214	txt=""
215	write(txt,"(a)") tnom(iq)
216	write(,)" tracer:",trim(txt)
217	! CO2
218	if (txt.eq."co2") then
219	q(:,iq)=0.95 ! kg /kg of atmosphere
220	qsurf(iq)=0. ! kg/m2 (not used for CO2)
221	! even better, look for a "profile_co2" input file
222	open(91,file='profile_co2',status='old',
223	& form='formatted',iostat=ierr)
224	if (ierr.eq.0) then
225	read(91,*) qsurf(iq)
226	do ilayer=1,nlayermx
227	read(91,*) q(ilayer,iq)
228	enddo
229	endif
230	close(91)
231	endif ! of if (txt.eq."co2")
232	! Allow for an initial profile of argon
233	! Can also be used to introduce a decaying tracer
234	! in the 1D (TBD) to study thermals
235	if (txt.eq."ar") then
236	!look for a "profile_ar" input file
237	open(91,file='profile_ar',status='old',
238	& form='formatted',iostat=ierr)
239	if (ierr.eq.0) then
240	read(91,*) qsurf(iq)
241	do ilayer=1,nlayermx
242	read(91,*) q(ilayer,iq)
243	enddo
244	else
245	write(,) "No profile_ar file!"
246	endif
247	close(91)
248	endif ! of if (txt.eq."ar")
249
250	! WATER VAPOUR
251	if (txt.eq."h2o_vap") then
252	!look for a "profile_h2o_vap" input file
253	open(91,file='profile_h2o_vap',status='old',
254	& form='formatted',iostat=ierr)
255	if (ierr.eq.0) then
256	read(91,*) qsurf(iq)
257	do ilayer=1,nlayermx
258	read(91,*) q(ilayer,iq)
259	enddo
260	else
261	write(,) "No profile_h2o_vap file!"
262	endif
263	close(91)
264	endif ! of if (txt.eq."h2o_ice")
265	! WATER ICE
266	if (txt.eq."h2o_ice") then
267	!look for a "profile_h2o_vap" input file
268	open(91,file='profile_h2o_ice',status='old',
269	& form='formatted',iostat=ierr)
270	if (ierr.eq.0) then
271	read(91,*) qsurf(iq)
272	do ilayer=1,nlayermx
273	read(91,*) q(ilayer,iq)
274	enddo
275	else
276	write(,) "No profile_h2o_ice file!"
277	endif
278	close(91)
279	endif ! of if (txt.eq."h2o_ice")
280	! DUST
281	!if (txt(1:4).eq."dust") then
282	! q(:,iq)=0.4 ! kg/kg of atmosphere
283	! qsurf(iq)=100 ! kg/m2
284	!endif
285	! DUST MMR
286	if (txt.eq."dust_mass") then
287	!look for a "profile_dust_mass" input file
288	open(91,file='profile_dust_mass',status='old',
289	& form='formatted',iostat=ierr)
290	if (ierr.eq.0) then
291	read(91,*) qsurf(iq)
292	do ilayer=1,nlayermx
293	read(91,*) q(ilayer,iq)
294	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
295	enddo
296	else
297	write(,) "No profile_dust_mass file!"
298	endif
299	close(91)
300	endif ! of if (txt.eq."dust_mass")
301	! DUST NUMBER
302	if (txt.eq."dust_number") then
303	!look for a "profile_dust_number" input file
304	open(91,file='profile_dust_number',status='old',
305	& form='formatted',iostat=ierr)
306	if (ierr.eq.0) then
307	read(91,*) qsurf(iq)
308	do ilayer=1,nlayermx
309	read(91,*) q(ilayer,iq)
310	enddo
311	else
312	write(,) "No profile_dust_number file!"
313	endif
314	close(91)
315	endif ! of if (txt.eq."dust_number")
316	! NB: some more initializations (chemistry) is done later
317	! CCN MASS
318	if (txt.eq."ccn_mass") then
319	!look for a "profile_ccn_mass" input file
320	open(91,file='profile_ccn_mass',status='old',
321	& form='formatted',iostat=ierr)
322	if (ierr.eq.0) then
323	read(91,*) qsurf(iq)
324	do ilayer=1,nlayermx
325	read(91,*) q(ilayer,iq)
326	enddo
327	else
328	write(,) "No profile_ccn_mass file!"
329	endif
330	close(91)
331	endif ! of if (txt.eq."ccn_mass")
332	! CCN NUMBER
333	if (txt.eq."ccn_number") then
334	!look for a "profile_ccn_number" input file
335	open(91,file='profile_ccn_number',status='old',
336	& form='formatted',iostat=ierr)
337	if (ierr.eq.0) then
338	read(91,*) qsurf(iq)
339	do ilayer=1,nlayermx
340	read(91,*) q(ilayer,iq)
341	enddo
342	else
343	write(,) "No profile_ccn_number file!"
344	endif
345	close(91)
346	endif ! of if (txt.eq."ccn_number")
347	enddo ! of do iq=1,nqmx
348
349	else
350	! we still need to set (dummy) tracer names for physdem1
351	nq=nqmx
352	do iq=1,nq
353	write(str7,'(a1,i2.2)')'q',iq
354	tnom(iq)=str7
355	enddo
356	endif ! of if (tracer)
357
358	!write(,) "testphys1d q", q(1,:)
359	!write(,) "testphys1d qsurf", qsurf
360
361	c Date et heure locale du debut du run
362	c ------------------------------------
363	c Date (en sols depuis le solstice de printemps) du debut du run
364	day0 = 0 ! default value for day0
365	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
366	call getin("day0",day0)
367	day=float(day0)
368	write(,) " day0 = ",day0
369	c Heure de demarrage
370	time=0 ! default value for time
371	write(,)'Initial local time (in hours, between 0 and 24)?'
372	call getin("time",time)
373	write(,)" time = ",time
374	time=time/24.E+0 ! convert time (hours) to fraction of sol
375
376	c Discretisation (Definition de la grille et des pas de temps)
377	c --------------
378	c
379	nlayer=nlayermx
380	nlevel=nlayer+1
381	nsoil=nsoilmx
382
383	day_step=48 ! default value for day_step
384	PRINT *,'Number of time steps per sol ?'
385	call getin("day_step",day_step)
386	write(,) " day_step = ",day_step
387
388	ndt=10 ! default value for ndt
389	PRINT *,'Number of sols to run ?'
390	call getin("ndt",ndt)
391	write(,) " ndt = ",ndt
392
393	ndt=ndt*day_step
394	dtphys=daysec/day_step
395	c Pression de surface sur la planete
396	c ------------------------------------
397	c
398	psurf=610. ! default value for psurf
399	PRINT *,'Surface pressure (Pa) ?'
400	call getin("psurf",psurf)
401	write(,) " psurf = ",psurf
402	c Pression de reference
403	pa=20.
404	preff=610.
405
406	c Proprietes des poussiere aerosol
407	c --------------------------------
408	tauvis=0.2 ! default value for tauvis
409	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
410	call getin("tauvis",tauvis)
411	write(,) " tauvis = ",tauvis
412
413	c latitude/longitude
414	c ------------------
415	lati(1)=0 ! default value for lati(1)
416	PRINT *,'latitude (in degrees) ?'
417	call getin("latitude",lati(1))
418	write(,) " latitude = ",lati(1)
419	lati(1)=lati(1)*pi/180.E+0
420	long(1)=0.E+0
421	long(1)=long(1)*pi/180.E+0
422
423	c Initialisation albedo / inertie du sol
424	c --------------------------------------
425	c
426	albedodat(1)=0.2 ! default value for albedodat
427	PRINT *,'Albedo of bare ground ?'
428	call getin("albedo",albedodat(1))
429	write(,) " albedo = ",albedodat(1)
430
431	inertiedat(1,1)=400 ! default value for inertiedat
432	PRINT *,'Soil thermal inertia (SI) ?'
433	call getin("inertia",inertiedat(1,1))
434	write(,) " inertia = ",inertiedat(1,1)
435
436	z0(1)=z0_default ! default value for roughness
437	write(,) 'Surface roughness length z0 (m)?'
438	call getin("z0",z0(1))
439	write(,) " z0 = ",z0(1)
440	c
441	c pour le schema d'ondes de gravite
442	c ---------------------------------
443	c
444	zmea(1)=0.E+0
445	zstd(1)=0.E+0
446	zsig(1)=0.E+0
447	zgam(1)=0.E+0
448	zthe(1)=0.E+0
449
450
451
452	c Initialisation speciales "physiq"
453	c ---------------------------------
454	c la surface de chaque maille est inutile en 1D --->
455	area(1)=1.E+0
456
457	c le geopotentiel au sol est inutile en 1D car tout est controle
458	c par la pression de surface --->
459	phisfi(1)=0.E+0
460
461	c "inifis" reproduit un certain nombre d'initialisations deja faites
462	c + lecture des clefs de callphys.def
463	c + calcul de la frequence d'appel au rayonnement
464	c + calcul des sinus et cosinus des longitude latitude
465
466	!Mars possible matter with dtphys in input and include!!!
467	CALL inifis(1,llm,day0,daysec,dtphys,
468	. lati,long,area,rad,g,r,cpp)
469	c Initialisation pour prendre en compte les vents en 1-D
470	c ------------------------------------------------------
471	ptif=2.E+0omegsinlat(1)
472
473	c vent geostrophique
474	gru=10. ! default value for gru
475	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
476	call getin("u",gru)
477	write(,) " u = ",gru
478	grv=0. !default value for grv
479	PRINT *,'meridional northward component of the geostrophic',
480	&' wind (m/s) ?'
481	call getin("v",grv)
482	write(,) " v = ",grv
483
484	c Initialisation des vents au premier pas de temps
485	DO ilayer=1,nlayer
486	u(ilayer)=gru
487	v(ilayer)=grv
488	ENDDO
489
490	c energie cinetique turbulente
491	DO ilevel=1,nlevel
492	q2(ilevel)=0.E+0
493	ENDDO
494
495	c glace de CO2 au sol
496	c -------------------
497	co2ice=0.E+0 ! default value for co2ice
498	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
499	call getin("co2ice",co2ice)
500	write(,) " co2ice = ",co2ice
501
502	c
503	c emissivite
504	c ----------
505	emis=emissiv
506	IF (co2ice.eq.1.E+0) THEN
507	emis=emisice(1) ! northern hemisphere
508	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
509	ENDIF
510
511
512
513	c calcul des pressions et altitudes en utilisant les niveaux sigma
514	c ----------------------------------------------------------------
515
516	c Vertical Coordinates
517	c """"""""""""""""""""
518	hybrid=.true.
519	PRINT *,'Hybrid coordinates ?'
520	call getin("hybrid",hybrid)
521	write(,) " hybrid = ", hybrid
522
523	CALL disvert
524
525	DO ilevel=1,nlevel
526	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
527	ENDDO
528
529	DO ilayer=1,nlayer
530	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
531	ENDDO
532
533	DO ilayer=1,nlayer
534	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
535	& /g
536	ENDDO
537
538
539	c profil de temperature au premier appel
540	c --------------------------------------
541	pks=psurf**rcp
542
543	c altitude en km dans profile: on divise zlay par 1000
544	tmp1(0)=0.E+0
545	DO ilayer=1,nlayer
546	tmp1(ilayer)=zlay(ilayer)/1000.E+0
547	ENDDO
548
549	call profile(nlayer+1,tmp1,tmp2)
550
551	tsurf=tmp2(0)
552	DO ilayer=1,nlayer
553	temp(ilayer)=tmp2(ilayer)
554	ENDDO
555
556
557
558	! Initialize soil properties and temperature
559	! ------------------------------------------
560	volcapa=1.e6 ! volumetric heat capacity
561	DO isoil=1,nsoil
562	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
563	tsoil(isoil)=tsurf ! soil temperature
564	ENDDO
565
566	! Initialize depths
567	! -----------------
568	do isoil=0,nsoil-1
569	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
570	enddo
571	do isoil=1,nsoil
572	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
573	enddo
574
575	c Initialisation des traceurs
576	c ---------------------------
577
578	if (photochem.or.callthermos) then
579	write(,) 'Especes chimiques initialisees'
580	! thermo=0: initialisation pour toutes les couches
581	thermo=0
582	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
583	$ thermo,qsurf)
584	endif
585
586	c Regarder si le sol est un reservoir de glace d'eau
587	c --------------------------------------------------
588	watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
589	print *,'Water ice cap on ground ?'
590	call getin("watercaptag",watercaptag)
591	write(,) " watercaptag = ",watercaptag
592
593
594	c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
595	c ----------------------------------------------------------------
596	c (effectuee avec "writeg1d" a partir de "physiq.F"
597	c ou tout autre subroutine physique, y compris celle ci).
598
599	g1d_nlayer=nlayer
600	g1d_nomfich='g1d.dat'
601	g1d_unitfich=40
602	g1d_nomctl='g1d.ctl'
603	g1d_unitctl=41
604	g1d_premier=.true.
605	g2d_premier=.true.
606
607	c Ecriture de "startfi"
608	c --------------------
609	c (Ce fichier sera aussitot relu au premier
610	c appel de "physiq", mais il est necessaire pour passer
611	c les variables purement physiques a "physiq"...
612
613	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
614	. dtphys,float(day0),time,tsurf,
615	. tsoil,co2ice,emis,q2,qsurf,area,albedodat,
616	. inertiedat,zmea,zstd,zsig,zgam,zthe)
617	c=======================================================================
618	c BOUCLE TEMPORELLE DU MODELE 1D
619	c=======================================================================
620	c
621	firstcall=.true.
622	lastcall=.false.
623
624	DO idt=1,ndt
625	c IF (idt.eq.ndt) lastcall=.true.
626	IF (idt.eq.ndt-day_step-1) then !test
627	lastcall=.true.
628	call solarlong(day*1.0,zls)
629	write(103,) 'Ls=',zls180./pi
630	write(103,) 'Lat=', lati(1)180./pi
631	write(103,) 'Tau=', tauvis/odprefpsurf
632	write(103,*) 'RunEnd - Atmos. Temp. File'
633	write(103,*) 'RunEnd - Atmos. Temp. File'
634	write(104,) 'Ls=',zls180./pi
635	write(104,*) 'Lat=', lati(1)
636	write(104,) 'Tau=', tauvis/odprefpsurf
637	write(104,*) 'RunEnd - Atmos. Temp. File'
638	ENDIF
639
640	c calcul du geopotentiel
641	c ~~~~~~~~~~~~~~~~~~~~~
642	DO ilayer=1,nlayer
643	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
644	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
645	ENDDO
646	phi(1)=pksh(1)(1.E+0-s(1))
647	DO ilayer=2,nlayer
648	phi(ilayer)=phi(ilayer-1)+
649	& pks(h(ilayer-1)+h(ilayer)).5E+0
650	& *(s(ilayer-1)-s(ilayer))
651
652	ENDDO
653
654	c appel de la physique
655	c --------------------
656	! write(,) "testphys1d avant q", q(1,:)
657	CALL physiq (1,llm,nqmx,
658	, firstcall,lastcall,
659	, day,time,dtphys,
660	, plev,play,phi,
661	, u, v,temp, q,
662	, w,
663	C - sorties
664	s du, dv, dtemp, dq,dpsurf,tracerdyn)
665	! write(,) "testphys1d apres q", q(1,:)
666
667
668	c evolution du vent : modele 1D
669	c -----------------------------
670
671	c la physique calcule les derivees temporelles de u et v.
672	c on y rajoute betement un effet Coriolis.
673	c
674	c DO ilayer=1,nlayer
675	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
676	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
677	c ENDDO
678
679	c Pour certain test : pas de coriolis a l'equateur
680	c if(lati(1).eq.0.) then
681	DO ilayer=1,nlayer
682	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
683	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
684	ENDDO
685	c end if
686	c
687	c
688	c Calcul du temps au pas de temps suivant
689	c ---------------------------------------
690	firstcall=.false.
691	time=time+dtphys/daysec
692	IF (time.gt.1.E+0) then
693	time=time-1.E+0
694	day=day+1
695	ENDIF
696
697	c calcul des vitesses et temperature au pas de temps suivant
698	c ----------------------------------------------------------
699
700	DO ilayer=1,nlayer
701	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
702	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
703	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
704	ENDDO
705
706	c calcul des pressions au pas de temps suivant
707	c ----------------------------------------------------------
708
709	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
710	DO ilevel=1,nlevel
711	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
712	ENDDO
713	DO ilayer=1,nlayer
714	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
715	ENDDO
716
717	! increment tracer values
718	DO iq = 1, nqmx
719	DO ilayer=1,nlayer
720	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
721	ENDDO
722	ENDDO
723
724	ENDDO ! fin de la boucle temporelle
725
726	c ========================================================
727	c GESTION DES SORTIE
728	c ========================================================
729
730	c fermeture pour conclure les sorties format grads dans "g1d.dat"
731	c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
732	c ou tout autre subroutine physique, y compris celle ci).
733
734	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
735	CALL endg1d(1,nlayer,zlay/1000.,ndt)
736
737	c ========================================================
738	END
739
740	c***********************************************************************
741	c***********************************************************************
742	c Subroutines Bidons utilise seulement en 3D, mais
743	c necessaire a la compilation de testphys1d en 1-D
744
745	subroutine gr_fi_dyn
746	RETURN
747	END
748
749	c***********************************************************************
750	c***********************************************************************
751
752	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format