Context Navigation

testphys1d.F @ 1621

Last change on this file since 1621 was 1621, checked in by emillour, 9 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added nqtot to tracer_h.F90.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

File size: 29.2 KB

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

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

Download in other formats:

Original Format