Context Navigation

testphys1d.F @ 1543

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

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

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

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

Download in other formats:

Original Format