Context Navigation

testphys1d.F @ 823

Last change on this file since 823 was 720, checked in by jbmadeleine, 12 years ago

Added the ability to set up the orbital parameters in run.def by

adding new calls to "getin" in testphys1d.F.

File size: 25.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 "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 Orbital parameters
414	c ------------------
415	print *,'Min. distance Sun-Mars (Mkm)?'
416	call getin("periheli",periheli)
417	write(,) " periheli = ",periheli
418
419	print *,'Max. distance Sun-Mars (Mkm)?'
420	call getin("aphelie",aphelie)
421	write(,) " aphelie = ",aphelie
422
423	print *,'Day of perihelion?'
424	call getin("periday",peri_day)
425	write(,) " periday = ",peri_day
426
427	print *,'Obliquity?'
428	call getin("obliquit",obliquit)
429	write(,) " obliquit = ",obliquit
430
431	c latitude/longitude
432	c ------------------
433	lati(1)=0 ! default value for lati(1)
434	PRINT *,'latitude (in degrees) ?'
435	call getin("latitude",lati(1))
436	write(,) " latitude = ",lati(1)
437	lati(1)=lati(1)*pi/180.E+0
438	long(1)=0.E+0
439	long(1)=long(1)*pi/180.E+0
440
441	c Initialisation albedo / inertie du sol
442	c --------------------------------------
443	c
444	albedodat(1)=0.2 ! default value for albedodat
445	PRINT *,'Albedo of bare ground ?'
446	call getin("albedo",albedodat(1))
447	write(,) " albedo = ",albedodat(1)
448
449	inertiedat(1,1)=400 ! default value for inertiedat
450	PRINT *,'Soil thermal inertia (SI) ?'
451	call getin("inertia",inertiedat(1,1))
452	write(,) " inertia = ",inertiedat(1,1)
453
454	z0(1)=z0_default ! default value for roughness
455	write(,) 'Surface roughness length z0 (m)?'
456	call getin("z0",z0(1))
457	write(,) " z0 = ",z0(1)
458	c
459	c pour le schema d'ondes de gravite
460	c ---------------------------------
461	c
462	zmea(1)=0.E+0
463	zstd(1)=0.E+0
464	zsig(1)=0.E+0
465	zgam(1)=0.E+0
466	zthe(1)=0.E+0
467
468
469
470	c Initialisation speciales "physiq"
471	c ---------------------------------
472	c la surface de chaque maille est inutile en 1D --->
473	area(1)=1.E+0
474
475	c le geopotentiel au sol est inutile en 1D car tout est controle
476	c par la pression de surface --->
477	phisfi(1)=0.E+0
478
479	c "inifis" reproduit un certain nombre d'initialisations deja faites
480	c + lecture des clefs de callphys.def
481	c + calcul de la frequence d'appel au rayonnement
482	c + calcul des sinus et cosinus des longitude latitude
483
484	!Mars possible matter with dtphys in input and include!!!
485	CALL inifis(1,llm,day0,daysec,dtphys,
486	. lati,long,area,rad,g,r,cpp)
487	c Initialisation pour prendre en compte les vents en 1-D
488	c ------------------------------------------------------
489	ptif=2.E+0omegsinlat(1)
490
491	c vent geostrophique
492	gru=10. ! default value for gru
493	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
494	call getin("u",gru)
495	write(,) " u = ",gru
496	grv=0. !default value for grv
497	PRINT *,'meridional northward component of the geostrophic',
498	&' wind (m/s) ?'
499	call getin("v",grv)
500	write(,) " v = ",grv
501
502	c Initialisation des vents au premier pas de temps
503	DO ilayer=1,nlayer
504	u(ilayer)=gru
505	v(ilayer)=grv
506	ENDDO
507
508	c energie cinetique turbulente
509	DO ilevel=1,nlevel
510	q2(ilevel)=0.E+0
511	ENDDO
512
513	c glace de CO2 au sol
514	c -------------------
515	co2ice=0.E+0 ! default value for co2ice
516	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
517	call getin("co2ice",co2ice)
518	write(,) " co2ice = ",co2ice
519
520	c
521	c emissivite
522	c ----------
523	emis=emissiv
524	IF (co2ice.eq.1.E+0) THEN
525	emis=emisice(1) ! northern hemisphere
526	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
527	ENDIF
528
529
530
531	c calcul des pressions et altitudes en utilisant les niveaux sigma
532	c ----------------------------------------------------------------
533
534	c Vertical Coordinates
535	c """"""""""""""""""""
536	hybrid=.true.
537	PRINT *,'Hybrid coordinates ?'
538	call getin("hybrid",hybrid)
539	write(,) " hybrid = ", hybrid
540
541	CALL disvert
542
543	DO ilevel=1,nlevel
544	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
545	ENDDO
546
547	DO ilayer=1,nlayer
548	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
549	ENDDO
550
551	DO ilayer=1,nlayer
552	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
553	& /g
554	ENDDO
555
556
557	c profil de temperature au premier appel
558	c --------------------------------------
559	pks=psurf**rcp
560
561	c altitude en km dans profile: on divise zlay par 1000
562	tmp1(0)=0.E+0
563	DO ilayer=1,nlayer
564	tmp1(ilayer)=zlay(ilayer)/1000.E+0
565	ENDDO
566
567	call profile(nlayer+1,tmp1,tmp2)
568
569	tsurf=tmp2(0)
570	DO ilayer=1,nlayer
571	temp(ilayer)=tmp2(ilayer)
572	ENDDO
573
574
575
576	! Initialize soil properties and temperature
577	! ------------------------------------------
578	volcapa=1.e6 ! volumetric heat capacity
579	DO isoil=1,nsoil
580	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
581	tsoil(isoil)=tsurf ! soil temperature
582	ENDDO
583
584	! Initialize depths
585	! -----------------
586	do isoil=0,nsoil-1
587	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
588	enddo
589	do isoil=1,nsoil
590	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
591	enddo
592
593	c Initialisation des traceurs
594	c ---------------------------
595
596	if (photochem.or.callthermos) then
597	write(,) 'Especes chimiques initialisees'
598	! thermo=0: initialisation pour toutes les couches
599	thermo=0
600	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
601	$ thermo,qsurf)
602	endif
603
604	c Regarder si le sol est un reservoir de glace d'eau
605	c --------------------------------------------------
606	watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
607	print *,'Water ice cap on ground ?'
608	call getin("watercaptag",watercaptag)
609	write(,) " watercaptag = ",watercaptag
610
611
612	c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
613	c ----------------------------------------------------------------
614	c (effectuee avec "writeg1d" a partir de "physiq.F"
615	c ou tout autre subroutine physique, y compris celle ci).
616
617	g1d_nlayer=nlayer
618	g1d_nomfich='g1d.dat'
619	g1d_unitfich=40
620	g1d_nomctl='g1d.ctl'
621	g1d_unitctl=41
622	g1d_premier=.true.
623	g2d_premier=.true.
624
625	c Ecriture de "startfi"
626	c --------------------
627	c (Ce fichier sera aussitot relu au premier
628	c appel de "physiq", mais il est necessaire pour passer
629	c les variables purement physiques a "physiq"...
630
631	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
632	. dtphys,float(day0),time,tsurf,
633	. tsoil,co2ice,emis,q2,qsurf,area,albedodat,
634	. inertiedat,zmea,zstd,zsig,zgam,zthe)
635	c=======================================================================
636	c BOUCLE TEMPORELLE DU MODELE 1D
637	c=======================================================================
638	c
639	firstcall=.true.
640	lastcall=.false.
641
642	DO idt=1,ndt
643	c IF (idt.eq.ndt) lastcall=.true.
644	IF (idt.eq.ndt-day_step-1) then !test
645	lastcall=.true.
646	call solarlong(day*1.0,zls)
647	write(103,) 'Ls=',zls180./pi
648	write(103,) 'Lat=', lati(1)180./pi
649	write(103,) 'Tau=', tauvis/odprefpsurf
650	write(103,*) 'RunEnd - Atmos. Temp. File'
651	write(103,*) 'RunEnd - Atmos. Temp. File'
652	write(104,) 'Ls=',zls180./pi
653	write(104,*) 'Lat=', lati(1)
654	write(104,) 'Tau=', tauvis/odprefpsurf
655	write(104,*) 'RunEnd - Atmos. Temp. File'
656	ENDIF
657
658	c calcul du geopotentiel
659	c ~~~~~~~~~~~~~~~~~~~~~
660	DO ilayer=1,nlayer
661	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
662	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
663	ENDDO
664	phi(1)=pksh(1)(1.E+0-s(1))
665	DO ilayer=2,nlayer
666	phi(ilayer)=phi(ilayer-1)+
667	& pks(h(ilayer-1)+h(ilayer)).5E+0
668	& *(s(ilayer-1)-s(ilayer))
669
670	ENDDO
671
672	c appel de la physique
673	c --------------------
674	! write(,) "testphys1d avant q", q(1,:)
675	CALL physiq (1,llm,nqmx,
676	, firstcall,lastcall,
677	, day,time,dtphys,
678	, plev,play,phi,
679	, u, v,temp, q,
680	, w,
681	C - sorties
682	s du, dv, dtemp, dq,dpsurf,tracerdyn)
683	! write(,) "testphys1d apres q", q(1,:)
684
685
686	c evolution du vent : modele 1D
687	c -----------------------------
688
689	c la physique calcule les derivees temporelles de u et v.
690	c on y rajoute betement un effet Coriolis.
691	c
692	c DO ilayer=1,nlayer
693	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
694	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
695	c ENDDO
696
697	c Pour certain test : pas de coriolis a l'equateur
698	c if(lati(1).eq.0.) then
699	DO ilayer=1,nlayer
700	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
701	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
702	ENDDO
703	c end if
704	c
705	c
706	c Calcul du temps au pas de temps suivant
707	c ---------------------------------------
708	firstcall=.false.
709	time=time+dtphys/daysec
710	IF (time.gt.1.E+0) then
711	time=time-1.E+0
712	day=day+1
713	ENDIF
714
715	c calcul des vitesses et temperature au pas de temps suivant
716	c ----------------------------------------------------------
717
718	DO ilayer=1,nlayer
719	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
720	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
721	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
722	ENDDO
723
724	c calcul des pressions au pas de temps suivant
725	c ----------------------------------------------------------
726
727	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
728	DO ilevel=1,nlevel
729	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
730	ENDDO
731	DO ilayer=1,nlayer
732	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
733	ENDDO
734
735	! increment tracer values
736	DO iq = 1, nqmx
737	DO ilayer=1,nlayer
738	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
739	ENDDO
740	ENDDO
741
742	ENDDO ! fin de la boucle temporelle
743
744	c ========================================================
745	c GESTION DES SORTIE
746	c ========================================================
747
748	c fermeture pour conclure les sorties format grads dans "g1d.dat"
749	c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
750	c ou tout autre subroutine physique, y compris celle ci).
751
752	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
753	CALL endg1d(1,nlayer,zlay/1000.,ndt)
754
755	c ========================================================
756	END
757
758	c***********************************************************************
759	c***********************************************************************
760	c Subroutines Bidons utilise seulement en 3D, mais
761	c necessaire a la compilation de testphys1d en 1-D
762
763	subroutine gr_fi_dyn
764	RETURN
765	END
766
767	c***********************************************************************
768	c***********************************************************************
769
770	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format