Context Navigation

testphys1d.F @ 224

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

Mars GCM:

Implemented using 'z0' roughness length map (important: 'z0' reference

field is in datafile surface.nc, which has also been updated).

made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc).
adapted 'tabfi.F' to use 'z0_default'.
adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used.
modified 'physdem1.F' to write 'z0' field to restart.nc
adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'.
adapted 'dustdevil.F' to use 'z0_default'.
'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant?").
modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file.

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

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

Download in other formats:

Original Format