Context Navigation

testphys1d.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 26.2 KB

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

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

Download in other formats:

Original Format