Context Navigation

testphys1d.F @ 317

Last change on this file since 317 was 283, checked in by aslmd, 13 years ago

LMDZ.MARS:

AS: J'ai fait ce commit et fait a la main un merge avec les modifications entre temps

24/08/11 == TN

Attempts to tune the water cycle by adding outliers

+ A few structural changes !!

watercap.h is now obsolete and removed -- all is in surfdat.h
watercaptag initialized in surfini.F (up to 7 areas defined) instead of initracer.F
- settings proposed by AS commented
- experiments by TN decommented. use with caution.
water ice albedo and thermal inertia in callphys.def and inifis.F
water ice albedo in surfini.F
water ice albedo computation in albedocaps.F90
alb_surfice is now obsolete in physiq.F, albedo_h2o_ice is used instead
frost_albedo_threshold defined in surfdat.h
water ice thermal inertia in soil.F

TODO: * calibrate thermal inertia and ice albedo

have a look at subgrid-scale ice with dryness ?

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

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

Download in other formats:

Original Format