Context Navigation

testphys1d.F @ 899

Last change on this file since 899 was 899, checked in by emillour, 12 years ago

Mars GCM:

added possibility to set slope parameters in testphys1d

File size: 25.9 KB

Line
1
2	PROGRAM testphys1d
3	! to use 'getin'
4	USE ioipsl_getincom
5	IMPLICIT NONE
6
7	c=======================================================================
8	c subject:
9	c --------
10	c PROGRAM useful to run physical part of the martian GCM in a 1D column
11	c
12	c Can be compiled with a command like (e.g. for 25 layers)
13	c "makegcm -p mars -d 25 testphys1d"
14	c It requires the files "testphys1d.def" "callphys.def"
15	c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
16	c and a file describing the sigma layers (e.g. "z2sig.def")
17	c
18	c author: Frederic Hourdin, R.Fournier,F.Forget
19	c -------
20	c
21	c update: 12/06/2003 including chemistry (S. Lebonnois)
22	c and water ice (F. Montmessin)
23	c
24	c=======================================================================
25
26	#include "dimensions.h"
27	#include "dimphys.h"
28	#include "dimradmars.h"
29	#include "comgeomfi.h"
30	#include "surfdat.h"
31	#include "slope.h"
32	#include "comsoil.h"
33	#include "comdiurn.h"
34	#include "callkeys.h"
35	#include "comcstfi.h"
36	#include "planete.h"
37	#include "comsaison.h"
38	#include "yomaer.h"
39	#include "control.h"
40	#include "comvert.h"
41	#include "netcdf.inc"
42	#include "comg1d.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	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
56	c
57	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
58	REAL day ! date durant le run
59	REAL time ! time (0<time<1 ; time=0.5 a midi)
60	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
61	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
62	REAL psurf,tsurf
63	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
64	REAL gru,grv ! prescribed "geostrophic" background wind
65	REAL temp(nlayermx) ! temperature at the middle of the layers
66	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
67	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
68	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
69	REAL co2ice ! co2ice layer (kg.m-2)
70	REAL emis ! surface layer
71	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
72	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
73
74	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
75	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
76	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
77	REAL dpsurf
78	REAL dq(nlayermx,nqmx)
79	REAL dqdyn(nlayermx,nqmx)
80
81	c Various intermediate variables
82	INTEGER thermo
83	REAL zls
84	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
85	REAL pks, ptif, w(nlayermx)
86	REAL qtotinit, mqtot(nqmx),qtot
87	INTEGER ierr, aslun
88	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
89	Logical tracerdyn
90	integer :: nq=1 ! number of tracers
91
92	character*2 str2
93	character (len=7) :: str7
94	character(len=44) :: txt
95
96	c=======================================================================
97
98	c=======================================================================
99	c INITIALISATION
100	c=======================================================================
101
102	c ------------------------------------------------------
103	c Prescribed constants to be set here
104	c ------------------------------------------------------
105
106	pi=2.E+0*asin(1.E+0)
107
108	c Mars planetary constants
109	c ----------------------------
110	rad=3397200. ! mars radius (m) ~3397200 m
111	daysec=88775. ! length of a sol (s) ~88775 s
112	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
113	g=3.72 ! gravity (m.s-2) ~3.72
114	mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49
115	rcp=.256793 ! = r/cp ~0.256793
116	r= 8.314511E+0 *1000.E+0/mugaz
117	cpp= r/rcp
118	year_day = 669 ! lenght of year (sols) ~668.6
119	periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66
120	aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22
121	peri_day = 485. ! perihelion date (sols since N. Spring)
122	obliquit = 25.2 ! Obliquity (deg) ~25.2
123
124	c Planetary Boundary Layer and Turbulence parameters
125	c --------------------------------------------------
126	z0_default = 1.e-2 ! surface roughness (m) ~0.01
127	emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8
128	lmixmin = 30 ! mixing length ~100
129
130	c cap properties and surface emissivities
131	c ----------------------------------------------------
132	emissiv= 0.95 ! Bare ground emissivity ~.95
133	emisice(1)=0.95 ! Northern cap emissivity
134	emisice(2)=0.95 ! Southern cap emisssivity
135	albedice(1)=0.5 ! Northern cap albedo
136	albedice(2)=0.5 ! Southern cap albedo
137	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
138	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
139	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
140	dtemisice(2) = 2. ! time scale for snow metamorphism (south
141
142
143	c ------------------------------------------------------
144	c Loading run parameters from "run.def" file
145	c ------------------------------------------------------
146
147
148	! check if 'run.def' file is around (otherwise reading parameters
149	! from callphys.def via getin() routine won't work.
150	open(99,file='run.def',status='old',form='formatted',
151	& iostat=ierr)
152	if (ierr.ne.0) then
153	write(,) 'Cannot find required file "run.def"'
154	write(,) ' (which should contain some input parameters'
155	write(,) ' along with the following line:'
156	write(,) ' INCLUDEDEF=callphys.def'
157	write(,) ' )'
158	write(,) ' ... might as well stop here ...'
159	stop
160	else
161	close(99)
162	endif
163
164	! check if we are going to run with or without tracers
165	write(,) "Run with or without tracer transport ?"
166	tracer=.false. ! default value
167	call getin("tracer",tracer)
168	write(,) " tracer = ",tracer
169
170	! while we're at it, check if there is a 'traceur.def' file
171	! and preocess it, if necessary. Otherwise initialize tracer names
172	if (tracer) then
173	! load tracer names from file 'traceur.def'
174	open(90,file='traceur.def',status='old',form='formatted',
175	& iostat=ierr)
176	if (ierr.ne.0) then
177	write(,) 'Cannot find required file "traceur.def"'
178	write(,) ' If you want to run with tracers, I need it'
179	write(,) ' ... might as well stop here ...'
180	stop
181	else
182	write(,) "testphys1d: Reading file traceur.def"
183	! read number of tracers:
184	read(90,*,iostat=ierr) nq
185	if (ierr.ne.0) then
186	write(,) "testphys1d: error reading number of tracers"
187	write(,) " (first line of traceur.def) "
188	stop
189	else
190	! check that the number of tracers is indeed nqmx
191	if (nq.ne.nqmx) then
192	write(,) "testphys1d: error, wrong number of tracers:"
193	write(,) "nq=",nq," whereas nqmx=",nqmx
194	stop
195	endif
196	endif
197	endif
198	! read tracer names from file traceur.def
199	do iq=1,nqmx
200	read(90,*,iostat=ierr) tnom(iq)
201	if (ierr.ne.0) then
202	write(,) 'testphys1d: error reading tracer names...'
203	stop
204	endif
205	enddo
206	close(90)
207
208	! initialize tracers here:
209	write(,) "testphys1d: initializing tracers"
210	q(:,:)=0 ! default, set everything to zero
211	qsurf(:)=0
212	! "smarter" initialization of some tracers
213	! (get values from "profile_*" files, if these are available)
214	do iq=1,nqmx
215	txt=""
216	write(txt,"(a)") tnom(iq)
217	write(,)" tracer:",trim(txt)
218	! CO2
219	if (txt.eq."co2") then
220	q(:,iq)=0.95 ! kg /kg of atmosphere
221	qsurf(iq)=0. ! kg/m2 (not used for CO2)
222	! even better, look for a "profile_co2" input file
223	open(91,file='profile_co2',status='old',
224	& form='formatted',iostat=ierr)
225	if (ierr.eq.0) then
226	read(91,*) qsurf(iq)
227	do ilayer=1,nlayermx
228	read(91,*) q(ilayer,iq)
229	enddo
230	endif
231	close(91)
232	endif ! of if (txt.eq."co2")
233	! Allow for an initial profile of argon
234	! Can also be used to introduce a decaying tracer
235	! in the 1D (TBD) to study thermals
236	if (txt.eq."ar") then
237	!look for a "profile_ar" input file
238	open(91,file='profile_ar',status='old',
239	& form='formatted',iostat=ierr)
240	if (ierr.eq.0) then
241	read(91,*) qsurf(iq)
242	do ilayer=1,nlayermx
243	read(91,*) q(ilayer,iq)
244	enddo
245	else
246	write(,) "No profile_ar file!"
247	endif
248	close(91)
249	endif ! of if (txt.eq."ar")
250
251	! WATER VAPOUR
252	if (txt.eq."h2o_vap") then
253	!look for a "profile_h2o_vap" input file
254	open(91,file='profile_h2o_vap',status='old',
255	& form='formatted',iostat=ierr)
256	if (ierr.eq.0) then
257	read(91,*) qsurf(iq)
258	do ilayer=1,nlayermx
259	read(91,*) q(ilayer,iq)
260	enddo
261	else
262	write(,) "No profile_h2o_vap file!"
263	endif
264	close(91)
265	endif ! of if (txt.eq."h2o_ice")
266	! WATER ICE
267	if (txt.eq."h2o_ice") then
268	!look for a "profile_h2o_vap" input file
269	open(91,file='profile_h2o_ice',status='old',
270	& form='formatted',iostat=ierr)
271	if (ierr.eq.0) then
272	read(91,*) qsurf(iq)
273	do ilayer=1,nlayermx
274	read(91,*) q(ilayer,iq)
275	enddo
276	else
277	write(,) "No profile_h2o_ice file!"
278	endif
279	close(91)
280	endif ! of if (txt.eq."h2o_ice")
281	! DUST
282	!if (txt(1:4).eq."dust") then
283	! q(:,iq)=0.4 ! kg/kg of atmosphere
284	! qsurf(iq)=100 ! kg/m2
285	!endif
286	! DUST MMR
287	if (txt.eq."dust_mass") then
288	!look for a "profile_dust_mass" input file
289	open(91,file='profile_dust_mass',status='old',
290	& form='formatted',iostat=ierr)
291	if (ierr.eq.0) then
292	read(91,*) qsurf(iq)
293	do ilayer=1,nlayermx
294	read(91,*) q(ilayer,iq)
295	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
296	enddo
297	else
298	write(,) "No profile_dust_mass file!"
299	endif
300	close(91)
301	endif ! of if (txt.eq."dust_mass")
302	! DUST NUMBER
303	if (txt.eq."dust_number") then
304	!look for a "profile_dust_number" input file
305	open(91,file='profile_dust_number',status='old',
306	& form='formatted',iostat=ierr)
307	if (ierr.eq.0) then
308	read(91,*) qsurf(iq)
309	do ilayer=1,nlayermx
310	read(91,*) q(ilayer,iq)
311	enddo
312	else
313	write(,) "No profile_dust_number file!"
314	endif
315	close(91)
316	endif ! of if (txt.eq."dust_number")
317	! NB: some more initializations (chemistry) is done later
318	! CCN MASS
319	if (txt.eq."ccn_mass") then
320	!look for a "profile_ccn_mass" input file
321	open(91,file='profile_ccn_mass',status='old',
322	& form='formatted',iostat=ierr)
323	if (ierr.eq.0) then
324	read(91,*) qsurf(iq)
325	do ilayer=1,nlayermx
326	read(91,*) q(ilayer,iq)
327	enddo
328	else
329	write(,) "No profile_ccn_mass file!"
330	endif
331	close(91)
332	endif ! of if (txt.eq."ccn_mass")
333	! CCN NUMBER
334	if (txt.eq."ccn_number") then
335	!look for a "profile_ccn_number" input file
336	open(91,file='profile_ccn_number',status='old',
337	& form='formatted',iostat=ierr)
338	if (ierr.eq.0) then
339	read(91,*) qsurf(iq)
340	do ilayer=1,nlayermx
341	read(91,*) q(ilayer,iq)
342	enddo
343	else
344	write(,) "No profile_ccn_number file!"
345	endif
346	close(91)
347	endif ! of if (txt.eq."ccn_number")
348	enddo ! of do iq=1,nqmx
349
350	else
351	! we still need to set (dummy) tracer names for physdem1
352	nq=nqmx
353	do iq=1,nq
354	write(str7,'(a1,i2.2)')'q',iq
355	tnom(iq)=str7
356	enddo
357	! and just to be clean, also initialize tracers to zero for physdem1
358	q(:,:)=0
359	qsurf(:)=0
360	endif ! of if (tracer)
361
362	!write(,) "testphys1d q", q(1,:)
363	!write(,) "testphys1d qsurf", qsurf
364
365	c Date and local time at beginning of run
366	c ---------------------------------------
367	c Date (in sols since spring solstice) at beginning of run
368	day0 = 0 ! default value for day0
369	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
370	call getin("day0",day0)
371	day=float(day0)
372	write(,) " day0 = ",day0
373	c Local time at beginning of run
374	time=0 ! default value for time
375	write(,)'Initial local time (in hours, between 0 and 24)?'
376	call getin("time",time)
377	write(,)" time = ",time
378	time=time/24.E+0 ! convert time (hours) to fraction of sol
379
380	c Discretization (Definition of grid and time steps)
381	c --------------
382	c
383	nlayer=nlayermx
384	nlevel=nlayer+1
385	nsoil=nsoilmx
386
387	day_step=48 ! default value for day_step
388	PRINT *,'Number of time steps per sol ?'
389	call getin("day_step",day_step)
390	write(,) " day_step = ",day_step
391
392	ndt=10 ! default value for ndt
393	PRINT *,'Number of sols to run ?'
394	call getin("ndt",ndt)
395	write(,) " ndt = ",ndt
396
397	ndt=ndt*day_step
398	dtphys=daysec/day_step
399
400	c Imposed surface pressure
401	c ------------------------------------
402	c
403	psurf=610. ! default value for psurf
404	PRINT *,'Surface pressure (Pa) ?'
405	call getin("psurf",psurf)
406	write(,) " psurf = ",psurf
407	c Reference pressures
408	pa=20. ! transition pressure (for hybrid coord.)
409	preff=610. ! reference surface pressure
410
411	c Aerosol properties
412	c --------------------------------
413	tauvis=0.2 ! default value for tauvis (dust opacity)
414	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
415	call getin("tauvis",tauvis)
416	write(,) " tauvis = ",tauvis
417
418	c Orbital parameters
419	c ------------------
420	print *,'Min. distance Sun-Mars (Mkm)?'
421	call getin("periheli",periheli)
422	write(,) " periheli = ",periheli
423
424	print *,'Max. distance Sun-Mars (Mkm)?'
425	call getin("aphelie",aphelie)
426	write(,) " aphelie = ",aphelie
427
428	print *,'Day of perihelion?'
429	call getin("periday",peri_day)
430	write(,) " periday = ",peri_day
431
432	print *,'Obliquity?'
433	call getin("obliquit",obliquit)
434	write(,) " obliquit = ",obliquit
435
436	c latitude/longitude
437	c ------------------
438	lati(1)=0 ! default value for lati(1)
439	PRINT *,'latitude (in degrees) ?'
440	call getin("latitude",lati(1))
441	write(,) " latitude = ",lati(1)
442	lati(1)=lati(1)*pi/180.E+0
443	long(1)=0.E+0
444	long(1)=long(1)*pi/180.E+0
445
446	c Initialize albedo / soil thermal inertia
447	c ----------------------------------------
448	c
449	albedodat(1)=0.2 ! default value for albedodat
450	PRINT *,'Albedo of bare ground ?'
451	call getin("albedo",albedodat(1))
452	write(,) " albedo = ",albedodat(1)
453
454	inertiedat(1,1)=400 ! default value for inertiedat
455	PRINT *,'Soil thermal inertia (SI) ?'
456	call getin("inertia",inertiedat(1,1))
457	write(,) " inertia = ",inertiedat(1,1)
458
459	z0(1)=z0_default ! default value for roughness
460	write(,) 'Surface roughness length z0 (m)?'
461	call getin("z0",z0(1))
462	write(,) " z0 = ",z0(1)
463
464	! Initialize local slope parameters (only matters if "callslope"
465	! is .true. in callphys.def)
466	! slope inclination angle (deg) 0: horizontal, 90: vertical
467	theta_sl(1)=0.0 ! default: no inclination
468	call getin("slope_inclination",theta_sl(1))
469	! slope orientation (deg)
470	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
471	psi_sl(1)=0.0 ! default value
472	call getin("slope_orientation",psi_sl(1))
473
474	c
475	c for the gravity wave scheme
476	c ---------------------------------
477	c
478	zmea(1)=0.E+0
479	zstd(1)=0.E+0
480	zsig(1)=0.E+0
481	zgam(1)=0.E+0
482	zthe(1)=0.E+0
483
484
485	c Specific initializations for "physiq"
486	c -------------------------------------
487	c mesh surface (not a very usefull quantity in 1D)
488	area(1)=1.E+0
489
490	c surface geopotential is not used (or useful) since in 1D
491	c everything is controled by surface pressure
492	phisfi(1)=0.E+0
493
494	c "inifis" does some initializations (some of which have already been
495	c done above!) and loads parameters set in callphys.def
496
497	!Mars possible matter with dtphys in input and include!!!
498	CALL inifis(1,llm,day0,daysec,dtphys,
499	. lati,long,area,rad,g,r,cpp)
500
501	c Initialization to take into account prescribed winds
502	c ------------------------------------------------------
503	ptif=2.E+0omegsinlat(1)
504
505	c geostrophic wind
506	gru=10. ! default value for gru
507	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
508	call getin("u",gru)
509	write(,) " u = ",gru
510	grv=0. !default value for grv
511	PRINT *,'meridional northward component of the geostrophic',
512	&' wind (m/s) ?'
513	call getin("v",grv)
514	write(,) " v = ",grv
515
516	c Initialize winds for first time step
517	DO ilayer=1,nlayer
518	u(ilayer)=gru
519	v(ilayer)=grv
520	ENDDO
521
522	c Initialize turbulente kinetic energy
523	DO ilevel=1,nlevel
524	q2(ilevel)=0.E+0
525	ENDDO
526
527	c CO2 ice on the surface
528	c -------------------
529	co2ice=0.E+0 ! default value for co2ice
530	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
531	call getin("co2ice",co2ice)
532	write(,) " co2ice = ",co2ice
533
534	c
535	c emissivity
536	c ----------
537	emis=emissiv
538	IF (co2ice.eq.1.E+0) THEN
539	emis=emisice(1) ! northern hemisphere
540	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
541	ENDIF
542
543
544
545	c Compute pressures and altitudes of atmospheric levels
546	c ----------------------------------------------------------------
547
548	c Vertical Coordinates
549	c """"""""""""""""""""
550	hybrid=.true.
551	PRINT *,'Hybrid coordinates ?'
552	call getin("hybrid",hybrid)
553	write(,) " hybrid = ", hybrid
554
555	CALL disvert
556
557	DO ilevel=1,nlevel
558	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
559	ENDDO
560
561	DO ilayer=1,nlayer
562	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
563	ENDDO
564
565	DO ilayer=1,nlayer
566	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
567	& /g
568	ENDDO
569
570
571	c Initialize temperature profile
572	c --------------------------------------
573	pks=psurf**rcp
574
575	c altitude in km in profile: divide zlay by 1000
576	tmp1(0)=0.E+0
577	DO ilayer=1,nlayer
578	tmp1(ilayer)=zlay(ilayer)/1000.E+0
579	ENDDO
580
581	call profile(nlayer+1,tmp1,tmp2)
582
583	tsurf=tmp2(0)
584	DO ilayer=1,nlayer
585	temp(ilayer)=tmp2(ilayer)
586	ENDDO
587
588
589
590	! Initialize soil properties and temperature
591	! ------------------------------------------
592	volcapa=1.e6 ! volumetric heat capacity
593	DO isoil=1,nsoil
594	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
595	tsoil(isoil)=tsurf ! soil temperature
596	ENDDO
597
598	! Initialize depths
599	! -----------------
600	do isoil=0,nsoil-1
601	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
602	enddo
603	do isoil=1,nsoil
604	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
605	enddo
606
607	c Initialize traceurs
608	c ---------------------------
609
610	if (photochem.or.callthermos) then
611	write(,) 'Initializing chemical species'
612	! thermo=0: initialize over all atmospheric layers
613	thermo=0
614	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
615	$ thermo,qsurf)
616	endif
617
618	c Check if the surface is a water ice reservoir
619	c --------------------------------------------------
620	watercaptag(ngridmx)=.false. ! Default: no water ice reservoir
621	print *,'Water ice cap on ground ?'
622	call getin("watercaptag",watercaptag)
623	write(,) " watercaptag = ",watercaptag
624
625
626	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
627	c ----------------------------------------------------------------
628	c (output done in "writeg1d", typically called by "physiq.F")
629
630	g1d_nlayer=nlayer
631	g1d_nomfich='g1d.dat'
632	g1d_unitfich=40
633	g1d_nomctl='g1d.ctl'
634	g1d_unitctl=41
635	g1d_premier=.true.
636	g2d_premier=.true.
637
638	c Write a "startfi" file
639	c --------------------
640	c This file will be read during the first call to "physiq".
641	c It is needed to transfert physics variables to "physiq"...
642
643	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
644	. dtphys,float(day0),time,tsurf,
645	. tsoil,co2ice,emis,q2,qsurf,area,albedodat,
646	. inertiedat,zmea,zstd,zsig,zgam,zthe)
647
648	c=======================================================================
649	c 1D MODEL TIME STEPPING LOOP
650	c=======================================================================
651	c
652	firstcall=.true.
653	lastcall=.false.
654
655	DO idt=1,ndt
656	c IF (idt.eq.ndt) lastcall=.true.
657	IF (idt.eq.ndt-day_step-1) then !test
658	lastcall=.true.
659	call solarlong(day*1.0,zls)
660	write(103,) 'Ls=',zls180./pi
661	write(103,) 'Lat=', lati(1)180./pi
662	write(103,) 'Tau=', tauvis/odprefpsurf
663	write(103,*) 'RunEnd - Atmos. Temp. File'
664	write(103,*) 'RunEnd - Atmos. Temp. File'
665	write(104,) 'Ls=',zls180./pi
666	write(104,*) 'Lat=', lati(1)
667	write(104,) 'Tau=', tauvis/odprefpsurf
668	write(104,*) 'RunEnd - Atmos. Temp. File'
669	ENDIF
670
671	c compute geopotential
672	c ~~~~~~~~~~~~~~~~~~~~~
673	DO ilayer=1,nlayer
674	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
675	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
676	ENDDO
677	phi(1)=pksh(1)(1.E+0-s(1))
678	DO ilayer=2,nlayer
679	phi(ilayer)=phi(ilayer-1)+
680	& pks(h(ilayer-1)+h(ilayer)).5E+0
681	& *(s(ilayer-1)-s(ilayer))
682
683	ENDDO
684
685	c call physics
686	c --------------------
687	! write(,) "testphys1d avant q", q(1,:)
688	CALL physiq (1,llm,nqmx,
689	, firstcall,lastcall,
690	, day,time,dtphys,
691	, plev,play,phi,
692	, u, v,temp, q,
693	, w,
694	C - outputs
695	s du, dv, dtemp, dq,dpsurf,tracerdyn)
696	! write(,) "testphys1d apres q", q(1,:)
697
698
699	c wind increment : specific for 1D
700	c --------------------------------
701
702	c The physics compute the tendencies on u and v,
703	c here we just add Coriolos effect
704	c
705	c DO ilayer=1,nlayer
706	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
707	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
708	c ENDDO
709
710	c For some tests : No coriolis force at equator
711	c if(lati(1).eq.0.) then
712	DO ilayer=1,nlayer
713	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
714	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
715	ENDDO
716	c end if
717	c
718	c
719	c Compute time for next time step
720	c ---------------------------------------
721	firstcall=.false.
722	time=time+dtphys/daysec
723	IF (time.gt.1.E+0) then
724	time=time-1.E+0
725	day=day+1
726	ENDIF
727
728	c compute winds and temperature for next time step
729	c ----------------------------------------------------------
730
731	DO ilayer=1,nlayer
732	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
733	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
734	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
735	ENDDO
736
737	c compute pressure for next time step
738	c ----------------------------------------------------------
739
740	psurf=psurf+dtphys*dpsurf ! surface pressure change
741	DO ilevel=1,nlevel
742	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
743	ENDDO
744	DO ilayer=1,nlayer
745	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
746	ENDDO
747
748	! increment tracers
749	DO iq = 1, nqmx
750	DO ilayer=1,nlayer
751	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
752	ENDDO
753	ENDDO
754
755	ENDDO ! of idt=1,ndt ! end of time stepping loop
756
757	c ========================================================
758	c OUTPUTS
759	c ========================================================
760
761	c finalize and close grads files "g1d.dat" and "g1d.ctl"
762
763	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
764	CALL endg1d(1,nlayer,zlay/1000.,ndt)
765
766	c ========================================================
767	END
768
769	c***********************************************************************
770	c***********************************************************************
771	c Dummy subroutines used only in 3D, but required to
772	c compile testphys1d (to cleanly use writediagfi)
773
774	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
775
776	IMPLICIT NONE
777
778	INTEGER im,jm,ngrid,nfield
779	REAL pdyn(im,jm,nfield)
780	REAL pfi(ngrid,nfield)
781
782	if (ngrid.ne.1) then
783	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
784	stop
785	endif
786
787	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
788
789	end
790
791	c***********************************************************************
792	c***********************************************************************
793
794	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format