Context Navigation

testphys1d.F @ 520

Last change on this file since 520 was 520, checked in by tnavarro, 13 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging

File size: 24.8 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	! NB: some more initializations (chemistry) is done later
315	! CCN MASS
316	if (txt.eq."ccn_mass") then
317	!look for a "profile_ccn_mass" input file
318	open(91,file='profile_ccn_mass',status='old',
319	& form='formatted',iostat=ierr)
320	if (ierr.eq.0) then
321	read(91,*) qsurf(iq)
322	do ilayer=1,nlayermx
323	read(91,*) q(ilayer,iq)
324	enddo
325	else
326	write(,) "No profile_ccn_mass file!"
327	endif
328	close(91)
329	endif ! of if (txt.eq."ccn_mass")
330	! CCN NUMBER
331	if (txt.eq."ccn_number") then
332	!look for a "profile_ccn_number" input file
333	open(91,file='profile_ccn_number',status='old',
334	& form='formatted',iostat=ierr)
335	if (ierr.eq.0) then
336	read(91,*) qsurf(iq)
337	do ilayer=1,nlayermx
338	read(91,*) q(ilayer,iq)
339	enddo
340	else
341	write(,) "No profile_ccn_number file!"
342	endif
343	close(91)
344	endif ! of if (txt.eq."ccn_number")
345	enddo ! of do iq=1,nqmx
346
347	else
348	! we still need to set (dummy) tracer names for physdem1
349	nq=nqmx
350	do iq=1,nq
351	write(str7,'(a1,i2.2)')'q',iq
352	tnom(iq)=str7
353	enddo
354	endif ! of if (tracer)
355
356	!write(,) "testphys1d q", q(1,:)
357	!write(,) "testphys1d qsurf", qsurf
358
359	c Date et heure locale du debut du run
360	c ------------------------------------
361	c Date (en sols depuis le solstice de printemps) du debut du run
362	day0 = 0 ! default value for day0
363	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
364	call getin("day0",day0)
365	day=float(day0)
366	write(,) " day0 = ",day0
367	c Heure de demarrage
368	time=0 ! default value for time
369	write(,)'Initial local time (in hours, between 0 and 24)?'
370	call getin("time",time)
371	write(,)" time = ",time
372	time=time/24.E+0 ! convert time (hours) to fraction of sol
373
374	c Discretisation (Definition de la grille et des pas de temps)
375	c --------------
376	c
377	nlayer=nlayermx
378	nlevel=nlayer+1
379	nsoil=nsoilmx
380
381	day_step=48 ! default value for day_step
382	PRINT *,'Number of time steps per sol ?'
383	call getin("day_step",day_step)
384	write(,) " day_step = ",day_step
385
386	ndt=10 ! default value for ndt
387	PRINT *,'Number of sols to run ?'
388	call getin("ndt",ndt)
389	write(,) " ndt = ",ndt
390
391	ndt=ndt*day_step
392	dtphys=daysec/day_step
393	c Pression de surface sur la planete
394	c ------------------------------------
395	c
396	psurf=610. ! default value for psurf
397	PRINT *,'Surface pressure (Pa) ?'
398	call getin("psurf",psurf)
399	write(,) " psurf = ",psurf
400	c Pression de reference
401	pa=20.
402	preff=610.
403
404	c Proprietes des poussiere aerosol
405	c --------------------------------
406	tauvis=0.2 ! default value for tauvis
407	print *,'Reference dust opacity at 700 Pa ?'
408	call getin("tauvis",tauvis)
409	write(,) " tauvis = ",tauvis
410
411	c latitude/longitude
412	c ------------------
413	lati(1)=0 ! default value for lati(1)
414	PRINT *,'latitude (in degrees) ?'
415	call getin("latitude",lati(1))
416	write(,) " latitude = ",lati(1)
417	lati(1)=lati(1)*pi/180.E+0
418	long(1)=0.E+0
419	long(1)=long(1)*pi/180.E+0
420
421	c Initialisation albedo / inertie du sol
422	c --------------------------------------
423	c
424	albedodat(1)=0.2 ! default value for albedodat
425	PRINT *,'Albedo of bare ground ?'
426	call getin("albedo",albedodat(1))
427	write(,) " albedo = ",albedodat(1)
428
429	inertiedat(1,1)=400 ! default value for inertiedat
430	PRINT *,'Soil thermal inertia (SI) ?'
431	call getin("inertia",inertiedat(1,1))
432	write(,) " inertia = ",inertiedat(1,1)
433
434	z0(1)=z0_default ! default value for roughness
435	write(,) 'Surface roughness length z0 (m)?'
436	call getin("z0",z0(1))
437	write(,) " z0 = ",z0(1)
438	c
439	c pour le schema d'ondes de gravite
440	c ---------------------------------
441	c
442	zmea(1)=0.E+0
443	zstd(1)=0.E+0
444	zsig(1)=0.E+0
445	zgam(1)=0.E+0
446	zthe(1)=0.E+0
447
448
449
450	c Initialisation speciales "physiq"
451	c ---------------------------------
452	c la surface de chaque maille est inutile en 1D --->
453	area(1)=1.E+0
454
455	c le geopotentiel au sol est inutile en 1D car tout est controle
456	c par la pression de surface --->
457	phisfi(1)=0.E+0
458
459	c "inifis" reproduit un certain nombre d'initialisations deja faites
460	c + lecture des clefs de callphys.def
461	c + calcul de la frequence d'appel au rayonnement
462	c + calcul des sinus et cosinus des longitude latitude
463
464	!Mars possible matter with dtphys in input and include!!!
465	CALL inifis(1,llm,day0,daysec,dtphys,
466	. lati,long,area,rad,g,r,cpp)
467	c Initialisation pour prendre en compte les vents en 1-D
468	c ------------------------------------------------------
469	ptif=2.E+0omegsinlat(1)
470
471	c vent geostrophique
472	gru=10. ! default value for gru
473	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
474	call getin("u",gru)
475	write(,) " u = ",gru
476	grv=0. !default value for grv
477	PRINT *,'meridional northward component of the geostrophic',
478	&' wind (m/s) ?'
479	call getin("v",grv)
480	write(,) " v = ",grv
481
482	c Initialisation des vents au premier pas de temps
483	DO ilayer=1,nlayer
484	u(ilayer)=gru
485	v(ilayer)=grv
486	ENDDO
487
488	c energie cinetique turbulente
489	DO ilevel=1,nlevel
490	q2(ilevel)=0.E+0
491	ENDDO
492
493	c glace de CO2 au sol
494	c -------------------
495	co2ice=0.E+0 ! default value for co2ice
496	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
497	call getin("co2ice",co2ice)
498	write(,) " co2ice = ",co2ice
499
500	c
501	c emissivite
502	c ----------
503	emis=emissiv
504	IF (co2ice.eq.1.E+0) THEN
505	emis=emisice(1) ! northern hemisphere
506	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
507	ENDIF
508
509
510
511	c calcul des pressions et altitudes en utilisant les niveaux sigma
512	c ----------------------------------------------------------------
513
514	c Vertical Coordinates
515	c """"""""""""""""""""
516	hybrid=.true.
517	PRINT *,'Hybrid coordinates ?'
518	call getin("hybrid",hybrid)
519	write(,) " hybrid = ", hybrid
520
521	CALL disvert
522
523	DO ilevel=1,nlevel
524	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
525	ENDDO
526
527	DO ilayer=1,nlayer
528	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
529	ENDDO
530
531	DO ilayer=1,nlayer
532	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
533	& /g
534	ENDDO
535
536
537	c profil de temperature au premier appel
538	c --------------------------------------
539	pks=psurf**rcp
540
541	c altitude en km dans profile: on divise zlay par 1000
542	tmp1(0)=0.E+0
543	DO ilayer=1,nlayer
544	tmp1(ilayer)=zlay(ilayer)/1000.E+0
545	ENDDO
546
547	call profile(nlayer+1,tmp1,tmp2)
548
549	tsurf=tmp2(0)
550	DO ilayer=1,nlayer
551	temp(ilayer)=tmp2(ilayer)
552	ENDDO
553
554
555
556	! Initialize soil properties and temperature
557	! ------------------------------------------
558	volcapa=1.e6 ! volumetric heat capacity
559	DO isoil=1,nsoil
560	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
561	tsoil(isoil)=tsurf ! soil temperature
562	ENDDO
563
564	! Initialize depths
565	! -----------------
566	do isoil=0,nsoil-1
567	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
568	enddo
569	do isoil=1,nsoil
570	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
571	enddo
572
573	c Initialisation des traceurs
574	c ---------------------------
575
576	if (photochem.or.callthermos) then
577	write(,) 'Especes chimiques initialisees'
578	! thermo=0: initialisation pour toutes les couches
579	thermo=0
580	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
581	$ thermo,qsurf)
582	endif
583
584	c Regarder si le sol est un reservoir de glace d'eau
585	c --------------------------------------------------
586	watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
587	print *,'Water ice cap on ground ?'
588	call getin("watercaptag",watercaptag)
589	write(,) " watercaptag = ",watercaptag
590
591
592	c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
593	c ----------------------------------------------------------------
594	c (effectuee avec "writeg1d" a partir de "physiq.F"
595	c ou tout autre subroutine physique, y compris celle ci).
596
597	g1d_nlayer=nlayer
598	g1d_nomfich='g1d.dat'
599	g1d_unitfich=40
600	g1d_nomctl='g1d.ctl'
601	g1d_unitctl=41
602	g1d_premier=.true.
603	g2d_premier=.true.
604
605	c Ecriture de "startfi"
606	c --------------------
607	c (Ce fichier sera aussitot relu au premier
608	c appel de "physiq", mais il est necessaire pour passer
609	c les variables purement physiques a "physiq"...
610
611	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
612	. dtphys,float(day0),time,tsurf,
613	. tsoil,co2ice,emis,q2,qsurf,area,albedodat,
614	. inertiedat,zmea,zstd,zsig,zgam,zthe)
615	c=======================================================================
616	c BOUCLE TEMPORELLE DU MODELE 1D
617	c=======================================================================
618	c
619	firstcall=.true.
620	lastcall=.false.
621
622	DO idt=1,ndt
623	c IF (idt.eq.ndt) lastcall=.true.
624	IF (idt.eq.ndt-day_step-1) then !test
625	lastcall=.true.
626	call solarlong(day*1.0,zls)
627	write(103,) 'Ls=',zls180./pi
628	write(103,) 'Lat=', lati(1)180./pi
629	write(103,) 'Tau=', tauvis/700psurf
630	write(103,*) 'RunEnd - Atmos. Temp. File'
631	write(103,*) 'RunEnd - Atmos. Temp. File'
632	write(104,) 'Ls=',zls180./pi
633	write(104,*) 'Lat=', lati(1)
634	write(104,) 'Tau=', tauvis/700psurf
635	write(104,*) 'RunEnd - Atmos. Temp. File'
636	ENDIF
637
638	c calcul du geopotentiel
639	c ~~~~~~~~~~~~~~~~~~~~~
640	DO ilayer=1,nlayer
641	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
642	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
643	ENDDO
644	phi(1)=pksh(1)(1.E+0-s(1))
645	DO ilayer=2,nlayer
646	phi(ilayer)=phi(ilayer-1)+
647	& pks(h(ilayer-1)+h(ilayer)).5E+0
648	& *(s(ilayer-1)-s(ilayer))
649
650	ENDDO
651
652	c appel de la physique
653	c --------------------
654	! write(,) "testphys1d avant q", q(1,:)
655	CALL physiq (1,llm,nqmx,
656	, firstcall,lastcall,
657	, day,time,dtphys,
658	, plev,play,phi,
659	, u, v,temp, q,
660	, w,
661	C - sorties
662	s du, dv, dtemp, dq,dpsurf,tracerdyn)
663	! write(,) "testphys1d apres q", q(1,:)
664
665
666	c evolution du vent : modele 1D
667	c -----------------------------
668
669	c la physique calcule les derivees temporelles de u et v.
670	c on y rajoute betement un effet Coriolis.
671	c
672	c DO ilayer=1,nlayer
673	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
674	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
675	c ENDDO
676
677	c Pour certain test : pas de coriolis a l'equateur
678	c if(lati(1).eq.0.) then
679	DO ilayer=1,nlayer
680	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
681	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
682	ENDDO
683	c end if
684	c
685	c
686	c Calcul du temps au pas de temps suivant
687	c ---------------------------------------
688	firstcall=.false.
689	time=time+dtphys/daysec
690	IF (time.gt.1.E+0) then
691	time=time-1.E+0
692	day=day+1
693	ENDIF
694
695	c calcul des vitesses et temperature au pas de temps suivant
696	c ----------------------------------------------------------
697
698	DO ilayer=1,nlayer
699	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
700	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
701	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
702	ENDDO
703
704	c calcul des pressions au pas de temps suivant
705	c ----------------------------------------------------------
706
707	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
708	DO ilevel=1,nlevel
709	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
710	ENDDO
711	DO ilayer=1,nlayer
712	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
713	ENDDO
714
715	! increment tracer values
716	DO iq = 1, nqmx
717	DO ilayer=1,nlayer
718	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
719	ENDDO
720	ENDDO
721
722	ENDDO ! fin de la boucle temporelle
723
724	c ========================================================
725	c GESTION DES SORTIE
726	c ========================================================
727
728	c fermeture pour conclure les sorties format grads dans "g1d.dat"
729	c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
730	c ou tout autre subroutine physique, y compris celle ci).
731
732	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
733	CALL endg1d(1,nlayer,zlay/1000.,ndt)
734
735	c ========================================================
736	END
737
738	c***********************************************************************
739	c***********************************************************************
740	c Subroutines Bidons utilise seulement en 3D, mais
741	c necessaire a la compilation de testphys1d en 1-D
742
743	subroutine gr_fi_dyn
744	RETURN
745	END
746
747	c***********************************************************************
748	c***********************************************************************
749
750	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format