Context Navigation

rcm1d.F @ 2236

Last change on this file since 2236 was 2116, checked in by jvatant, 6 years ago
Fix an error in tracer management in 1D when running with both mufi and chem. --JVO
Property svn:executable set to ``*
File size: 31.1 KB

Line
1	program rcm1d
2
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 tracer_h
9	use surfdat_h, only: albedodat, phisfi,
10	& zmea, zstd, zsig, zgam, zthe,
11	& emissiv, emisice, iceradius,
12	& dtemisice
13	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
14	use comchem_h, only: nkim, nlaykim_up, cnames, ykim_up
15	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
16	use phyredem, only: physdem0,physdem1
17	use geometry_mod, only: init_geometry
18	use planete_mod, only: apoastr,periastr,year_day,peri_day,
19	& obliquit,nres,z0,lmixmin,emin_turb,coefvis,coefir,
20	& timeperi,e_elips,p_elips
21	use comcstfi_mod, only: pi, cpp, rad, g, r,
22	& mugaz, rcp, omeg
23	use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy,
24	& nday, iphysiq
25	use callkeys_mod, only: tracer,check_cpp_match,rings_shadow,
26	& specOLR,pceil,callchim
27	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
28	& presnivs,pseudoalt,scaleheight
29	USE vertical_layers_mod, ONLY: init_vertical_layers
30	USE logic_mod, ONLY: hybrid,autozlevs
31	use regular_lonlat_mod, only: init_regular_lonlat
32	use planete_mod, only: ini_planete_mod
33	use physics_distribution_mod, only: init_physics_distribution
34	use regular_lonlat_mod, only: init_regular_lonlat
35	use mod_interface_dyn_phys, only: init_interface_dyn_phys
36	use inifis_mod, only: inifis
37	use physiq_mod, only: physiq
38	implicit none
39
40	!==================================================================
41	!
42	! Purpose
43	! -------
44	! Run the physics package of the universal model in a 1D column.
45	!
46	! It can be compiled with a command like (e.g. for 25 layers):
47	! "makegcm -p std -d 25 rcm1d"
48	! It requires the files "callphys.def", "z2sig.def",
49	! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
50	!
51	! Authors
52	! -------
53	! Frederic Hourdin
54	! R. Fournier
55	! F. Forget
56	! F. Montmessin (water ice added)
57	! R. Wordsworth
58	! B. Charnay
59	! A. Spiga
60	! J. Leconte (2012)
61	! J. Vatant d'Ollone (2017)
62	!
63	!==================================================================
64
65	#include "dimensions.h"
66	#include "paramet.h"
67	!include "dimphys.h"
68	#include "netcdf.inc"
69	#include "comgeom.h"
70
71	c --------------------------------------------------------------
72	c Declarations
73	c --------------------------------------------------------------
74	c
75	INTEGER unitstart ! unite d'ecriture de "startfi"
76	INTEGER nlayer,nlevel,nsoil,ndt
77	INTEGER ilayer,ilevel,isoil,idt,iq
78	LOGICAl firstcall,lastcall
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(llm) ! Pressure at the middle of the layers (Pa)
84	REAL plev(llm+1) ! intermediate pressure levels (pa)
85	REAL psurf,tsurf(1)
86	REAL u(llm),v(llm) ! zonal, meridional wind
87	REAL gru,grv ! prescribed "geostrophic" background wind
88	REAL temp(llm) ! 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,ALLOCATABLE :: tsoil(:) ! subsurface soil temperature (K)
92	REAL emis(1) ! surface layer
93	REAL q2(llm+1) ! Turbulent Kinetic Energy
94	REAL zlay(llm) ! altitude estimee dans les couches (km)
95
96	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
97	REAL du(llm),dv(llm),dtemp(llm)
98	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
99	REAL dpsurf(1)
100	REAL,ALLOCATABLE :: dq(:,:)
101	REAL,ALLOCATABLE :: dqdyn(:,:)
102
103	c Various intermediate variables
104	! INTEGER thermo
105	REAL zls
106	REAL phi(llm),h(llm),s(llm)
107	REAL pks, ptif, w(llm)
108	INTEGER ierr, aslun
109	REAL tmp1(0:llm),tmp2(0:llm)
110	integer :: nq !=1 ! number of tracers
111
112	character(len=44) :: txt
113
114	logical saveprofile
115
116	! added by RW for zlay computation
117	real Hscale, Hmax, rho, dz
118
119	! added by RW for autozlevs computation
120	real nu, xx, pMIN, zlev, Htop
121	real logplevs(llm)
122
123	! added by JVO
124	REAL tankCH4(1)
125
126	CHARACTER*38 file_prof
127	LOGICAL findprof
128
129	real :: latitude(1), longitude(1), cell_area(1)
130
131	c=======================================================================
132	c INITIALISATION
133	c=======================================================================
134	! initialize "serial/parallel" related stuff
135	! CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
136	! CALL init_phys_lmdz(1,1,llm,1,(/1/))
137	! call initcomgeomphy
138
139	!! those are defined in surfdat_h.F90
140	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(1))
141	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(1))
142	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(1))
143	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(1))
144	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(1))
145	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(1))
146	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(1))
147	!! those are defined in comdiurn_h.F90
148	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(1))
149	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(1))
150	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(1))
151	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(1))
152
153	saveprofile=.false.
154	saveprofile=.true.
155
156	c ----------------------------------------
157	c Default values (corresponding to Mars)
158	c ----------------------------------------
159
160	pi=2.E+0*asin(1.E+0)
161
162	c Parametres Couche limite et Turbulence
163	c --------------------------------------
164	z0 = 1.e-2 ! surface roughness (m) ~0.01
165	emin_turb = 1.e-6 ! energie minimale ~1.e-8
166	lmixmin = 30 ! longueur de melange ~100
167
168	c propriete optiques des calottes et emissivite du sol
169	c ----------------------------------------------------
170	emissiv= 0.95 ! Emissivite du sol martien ~.95
171	emisice(1)=0.95 ! Emissivite calotte nord
172	emisice(2)=0.95 ! Emissivite calotte sud
173
174	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
175	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
176	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
177	dtemisice(2) = 2. ! time scale for snow metamorphism (south
178	hybrid=.false.
179
180	c ------------------------------------------------------
181	c Load parameters from "run.def" and "gases.def"
182	c ------------------------------------------------------
183
184	! check if 'rcm1d.def' file is around
185	open(90,file='rcm1d.def',status='old',form='formatted',
186	& iostat=ierr)
187	if (ierr.ne.0) then
188	write(,) 'Cannot find required file "rcm1d.def"'
189	write(,) 'which should contain some input parameters'
190	write(,) ' ... might as well stop here ...'
191	stop
192	else
193	close(90)
194	endif
195
196	! now, run.def is needed anyway. so we create a dummy temporary one
197	! for ioipsl to work. if a run.def is already here, stop the
198	! program and ask the user to do a bit of cleaning
199	open(90,file='run.def',status='old',form='formatted',
200	& iostat=ierr)
201	if (ierr.eq.0) then
202	close(90)
203	write(,) 'There is already a run.def file.'
204	write(,) 'This is not compatible with 1D runs.'
205	write(,) 'Please remove the file and restart the run.'
206	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
207	stop
208	else
209	call system('touch run.def')
210	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
211	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
212	endif
213
214	! check if we are going to run with or without tracers
215	write(,) "Run with or without tracer transport ?"
216	tracer=.false. ! default value
217	call getin("tracer",tracer)
218	write(,) " tracer = ",tracer
219
220	! OK. now that run.def has been read once -- any variable is in memory.
221	! so we can dump the dummy run.def
222	! call system("rm -rf run.def") ! Ehouarn: delay this to after inifis
223
224	! while we're at it, check if there is a 'traceur.def' file
225	! and preocess it, if necessary. Otherwise initialize tracer names
226	if (tracer) then
227	! load tracer names from file 'traceur.def'
228	open(90,file='traceur.def',status='old',form='formatted',
229	& iostat=ierr)
230	if (ierr.eq.0) then
231	write(,) "rcm1d: Reading file traceur.def"
232	! read number of tracers:
233	read(90,*,iostat=ierr) nq
234	nqtot=nq ! set value of nqtot (in infotrac module) as nq
235	if (ierr.ne.0) then
236	write(,) "rcm1d: error reading number of tracers"
237	write(,) " (first line of traceur.def) "
238	stop
239	endif
240	if (nq>0) then
241	allocate(tname(nq))
242	allocate(q(llm,nq))
243	allocate(qsurf(nq))
244	allocate(dq(llm,nq))
245	allocate(dqdyn(llm,nq))
246	else
247	write(,) "rcm1d: Error. You set tracer=.true."
248	write(,) " but # of tracers in traceur.def is ",nq
249	stop
250	endif
251
252	do iq=1,nq
253	! minimal version, just read in the tracer names, 1 per line
254	read(90,*,iostat=ierr) tname(iq)
255	if (ierr.ne.0) then
256	write(,) 'rcm1d: error reading tracer names...'
257	stop
258	endif
259	enddo !of do iq=1,nq
260	else
261	write(,) 'Cannot find required file "traceur.def"'
262	write(,) ' If you want to run with tracers, I need it'
263	write(,) ' ... might as well stop here ...'
264	stop
265	endif
266	close(90)
267
268	else ! of if (tracer)
269	nqtot=0
270	nq=0
271	! still, make allocations for 1 dummy tracer
272	allocate(tname(1))
273	allocate(qsurf(1))
274	allocate(q(llm,1))
275	allocate(dq(llm,1))
276
277	! Check that tracer boolean is compliant with number of tracers
278	! -- otherwise there is an error (and more generally we have to be consistent)
279	if (nq .ge. 1) then
280	write(,) "------------------------------"
281	write(,) "rcm1d: You set tracer=.false."
282	write(,) " But set number of tracers to ",nq
283	write(,) " > If you want tracers, set tracer=.true."
284	write(,) "------------------------------"
285	stop
286	endif
287	endif ! of if (tracer)
288
289	!!! We have to check that check_cpp_match is F for 1D computations
290	!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
291	check_cpp_match = .false.
292	call getin("check_cpp_match",check_cpp_match)
293	if (check_cpp_match) then
294	print*,"In 1D modeling, check_cpp_match is supposed to be F"
295	print*,"Please correct callphys.def"
296	stop
297	endif
298
299	!!! GEOGRAPHICAL INITIALIZATIONS
300	!!! AREA. useless in 1D
301	cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
302	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
303	phisfi(1)=0.E+0
304	!!! LATITUDE. can be set.
305	latitude=0 ! default value for latitude
306	PRINT *,'latitude (in degrees) ?'
307	call getin("latitude",latitude)
308	write(,) " latitude = ",latitude
309	latitude=latitude*pi/180.E+0
310	!!! LONGITUDE. useless in 1D.
311	longitude=0.E+0
312	longitude=longitude*pi/180.E+0
313
314
315	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
316	!!!! PLANETARY CONSTANTS !!!!
317	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
318
319	g = -99999.
320	PRINT *,'GRAVITY in m s-2 ?'
321	call getin("g",g)
322	IF (g.eq.-99999.) THEN
323	PRINT *,"STOP. I NEED g IN RCM1D.DEF."
324	STOP
325	ELSE
326	PRINT *,"--> g = ",g
327	ENDIF
328
329	rad = -99999.
330	PRINT *,'PLANETARY RADIUS in m ?'
331	call getin("rad",rad)
332	! Planetary radius is needed to compute shadow of the rings
333	IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
334	PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
335	STOP
336	ELSE
337	PRINT *,"--> rad = ",rad
338	ENDIF
339
340	daysec = -99999.
341	PRINT *,'LENGTH OF A DAY in s ?'
342	call getin("daysec",daysec)
343	IF (daysec.eq.-99999.) THEN
344	PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
345	STOP
346	ELSE
347	PRINT *,"--> daysec = ",daysec
348	ENDIF
349	omeg=4.*asin(1.)/(daysec)
350	PRINT *,"OK. FROM THIS I WORKED OUT:"
351	PRINT *,"--> omeg = ",omeg
352
353	year_day = -99999.
354	PRINT *,'LENGTH OF A YEAR in days ?'
355	call getin("year_day",year_day)
356	IF (year_day.eq.-99999.) THEN
357	PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
358	STOP
359	ELSE
360	PRINT *,"--> year_day = ",year_day
361	ENDIF
362
363	periastr = -99999.
364	PRINT *,'MIN DIST STAR-PLANET in AU ?'
365	call getin("periastr",periastr)
366	IF (periastr.eq.-99999.) THEN
367	PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
368	STOP
369	ELSE
370	PRINT *,"--> periastr = ",periastr
371	ENDIF
372
373	apoastr = -99999.
374	PRINT *,'MAX DIST STAR-PLANET in AU ?'
375	call getin("apoastr",apoastr)
376	IF (apoastr.eq.-99999.) THEN
377	PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
378	STOP
379	ELSE
380	PRINT *,"--> apoastr = ",apoastr
381	ENDIF
382
383	peri_day = -99999.
384	PRINT *,'DATE OF PERIASTRON in days ?'
385	call getin("peri_day",peri_day)
386	IF (peri_day.eq.-99999.) THEN
387	PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
388	STOP
389	ELSE IF (peri_day.gt.year_day) THEN
390	PRINT *,"STOP. peri_day.gt.year_day"
391	STOP
392	ELSE
393	PRINT *,"--> peri_day = ", peri_day
394	ENDIF
395
396	obliquit = -99999.
397	PRINT *,'OBLIQUITY in deg ?'
398	call getin("obliquit",obliquit)
399	IF (obliquit.eq.-99999.) THEN
400	PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
401	STOP
402	ELSE
403	PRINT *,"--> obliquit = ",obliquit
404	ENDIF
405
406	psurf = -99999.
407	PRINT *,'SURFACE PRESSURE in Pa ?'
408	call getin("psurf",psurf)
409	IF (psurf.eq.-99999.) THEN
410	PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
411	STOP
412	ELSE
413	PRINT *,"--> psurf = ",psurf
414	ENDIF
415	!! we need reference pressures
416	pa=psurf/30.
417	preff=psurf ! these values are not needed in 1D (are you sure JL12)
418
419	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420	!!!! END PLANETARY CONSTANTS !!!!
421	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
422
423	c Date et heure locale du debut du run
424	c ------------------------------------
425	c Date (en sols depuis le solstice de printemps) du debut du run
426	day0 = 0 ! default value for day0
427	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
428	call getin("day0",day0)
429	day=float(day0)
430	write(,) " day0 = ",day0
431	c Heure de demarrage
432	time=0 ! default value for time
433	write(,)'Initial local time (in hours, between 0 and 24)?'
434	call getin("time",time)
435	write(,)" time = ",time
436	time=time/24.E+0 ! convert time (hours) to fraction of sol
437
438
439	c Discretisation (Definition de la grille et des pas de temps)
440	c --------------
441	c
442	nlayer=llm
443	nlevel=nlayer+1
444
445	day_step=48 ! default value for day_step
446	PRINT *,'Number of time steps per sol ?'
447	call getin("day_step",day_step)
448	write(,) " day_step = ",day_step
449
450	iphysiq=1 ! in 1D model physics are called evry time step
451	ecritphy=day_step ! default value for ecritphy
452	PRINT *,'Nunber of steps between writediagfi ?'
453	call getin("ecritphy",ecritphy)
454	write(,) " ecritphy = ",ecritphy
455
456	ndt=10 ! default value for ndt
457	PRINT *,'Number of sols to run ?'
458	call getin("ndt",ndt)
459	write(,) " ndt = ",ndt
460	nday=ndt
461
462	ndt=ndt*day_step
463	dtphys=daysec/day_step
464	write(,)"-------------------------------------"
465	write(,)"-------------------------------------"
466	write(,)"--> Physical timestep is ",dtphys
467	write(,)"-------------------------------------"
468	write(,)"-------------------------------------"
469
470	! initializations, as with iniphysiq.F90 for the 3D GCM
471	call init_physics_distribution(regular_lonlat,4,
472	& 1,1,1,nlayer,1)
473	call init_interface_dyn_phys
474	CALL init_regular_lonlat(1,1,longitude,latitude,
475	& (/0.,0./),(/0.,0./))
476	call init_geometry(1,longitude,latitude,
477	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
478	& cell_area)
479	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
480	! call init_vertical_layers(nlayer,preff,scaleheight,
481	! & ap,bp,aps,bps,presnivs,pseudoalt)
482	call init_dimphy(1,nlayer) ! Initialize dimphy module
483	call ini_planete_mod(nlayer,preff,ap,bp)
484
485	!!! CALL INIFIS
486	!!! - read callphys.def
487	!!! - calculate sine and cosine of longitude and latitude
488	!!! - calculate mugaz and cp
489	!!! NB: some operations are being done dummily in inifis in 1D
490	!!! - physical constants: nevermind, things are done allright below
491	!!! - physical frequency: nevermind, in inifis this is a simple print
492	cpp=-9999. ! dummy init for inifis, will be rewrite later on
493	r=-9999. ! dummy init for inifis, will be rewrite later on
494	CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
495	. latitude,longitude,cell_area,rad,g,r,cpp)
496
497	nsoil=nsoilmx
498	allocate(tsoil(nsoilmx))
499	!! those are defined in comsoil_h.F90
500	IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
501	IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
502	IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
503
504	! At this point, both getin() and getin_p() functions have been used,
505	! and the run.def file can be removed.
506	call system("rm -rf run.def")
507
508	!!! We check everything is OK.
509	PRINT *,"CHECK"
510	PRINT *,"--> mugaz = ",mugaz
511	PRINT *,"--> cpp = ",cpp
512	r = 8.314511E+0 * 1000.E+0 / mugaz
513	rcp = r / cpp
514
515	c output spectrum?
516	write(,)"Output spectral OLR?"
517	specOLR=.false.
518	call getin("specOLR",specOLR)
519	write(,)" specOLR = ",specOLR
520
521	c
522	c pour le schema d'ondes de gravite
523	c ---------------------------------
524	zmea(1)=0.E+0
525	zstd(1)=0.E+0
526	zsig(1)=0.E+0
527	zgam(1)=0.E+0
528	zthe(1)=0.E+0
529
530	c Initialisation des traceurs
531	c ---------------------------
532
533	DO iq=1,nq
534	DO ilayer=1,nlayer
535	q(ilayer,iq) = 0.
536	ENDDO
537	ENDDO
538
539	DO iq=1,nq
540	qsurf(iq) = 0.
541	ENDDO
542
543	if (tracer) then ! Initialize tracers here.
544
545	write(,) "rcm1d : initializing tracers profiles"
546
547	do iq=1,nq
548
549	txt=""
550	write(txt,"(a)") tname(iq)
551	write(,)" tracer:",trim(txt)
552
553	enddo ! of do iq=1,nq
554
555	endif ! of tracer
556
557	c Initialisation pour prendre en compte les vents en 1-D
558	c ------------------------------------------------------
559	ptif=2.E+0omegsinlat(1)
560
561
562	c vent geostrophique
563	gru=10. ! default value for gru
564	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
565	call getin("u",gru)
566	write(,) " u = ",gru
567	grv=0. !default value for grv
568	PRINT *,'meridional northward component of the geostrophic',
569	&' wind (m/s) ?'
570	call getin("v",grv)
571	write(,) " v = ",grv
572
573	! To be clean, also set vertical winds to zero
574	w(1:nlayer)=0
575
576	c Initialisation des vents au premier pas de temps
577	DO ilayer=1,nlayer
578	u(ilayer)=gru
579	v(ilayer)=grv
580	ENDDO
581
582	c energie cinetique turbulente
583	DO ilevel=1,nlevel
584	q2(ilevel)=0.E+0
585	ENDDO
586
587	c emissivity / surface co2 ice ( + h2o ice??)
588	c -------------------------------------------
589	emis(1)=emissiv ! default value for emissivity
590	PRINT *,'Emissivity of bare ground ?'
591	call getin("emis",emis(1))
592	write(,) " emis = ",emis(1)
593	emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
594	! value if there is no snow
595
596	c calcul des pressions et altitudes en utilisant les niveaux sigma
597	c ----------------------------------------------------------------
598
599	c Vertical Coordinates
600	c """"""""""""""""""""
601	hybrid=.true.
602	PRINT *,'Hybrid coordinates ?'
603	call getin("hybrid",hybrid)
604	write(,) " hybrid = ", hybrid
605
606
607	autozlevs=.false.
608	PRINT *,'Auto-discretise vertical levels ?'
609	call getin("autozlevs",autozlevs)
610	write(,) " autozlevs = ", autozlevs
611
612	pceil = psurf / 1000.0 ! Pascals
613	PRINT *,'Ceiling pressure (Pa) ?'
614	call getin("pceil",pceil)
615	write(,) " pceil = ", pceil
616
617	! Test of incompatibility:
618	! if autozlevs used, cannot have hybrid too
619
620	if (autozlevs.and.hybrid) then
621	print*,'Cannot use autozlevs and hybrid together!'
622	call abort
623	endif
624
625	if(autozlevs)then
626
627	open(91,file="z2sig.def",form='formatted')
628	read(91,*) Hscale
629	DO ilayer=1,nlayer-2
630	read(91,*) Hmax
631	enddo
632	close(91)
633
634
635	print*,'Hmax = ',Hmax,' km'
636	print*,'Auto-shifting Hscale to:'
637	! Hscale = Hmax / log(psurf/100.0)
638	Hscale = Hmax / log(psurf/pceil)
639	print*,'Hscale = ',Hscale,' km'
640
641	! none of this matters if we dont care about zlay
642
643	endif
644
645	call disvert
646	! now that disvert has been called, initialize module vertical_layers_mod
647	call init_vertical_layers(nlayer,preff,scaleheight,
648	& ap,bp,aps,bps,presnivs,pseudoalt)
649
650	if(.not.autozlevs)then
651	! we want only the scale height from z2sig, in order to compute zlay
652	open(91,file="z2sig.def",form='formatted')
653	read(91,*) Hscale
654	close(91)
655	endif
656
657	! if(autozlevs)then
658	! open(94,file="Hscale.temp",form='formatted')
659	! read(94,*) Hscale
660	! close(94)
661	! endif
662
663	DO ilevel=1,nlevel
664	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
665	ENDDO
666
667	DO ilayer=1,nlayer
668	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
669	ENDDO
670
671	DO ilayer=1,nlayer
672	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
673	! & /g
674	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
675	ENDDO
676
677	! endif
678
679	c profil de temperature au premier appel
680	c --------------------------------------
681	pks=psurf**rcp
682
683	c altitude en km dans profile: on divise zlay par 1000
684	tmp1(0)=0.E+0
685	DO ilayer=1,nlayer
686	tmp1(ilayer)=zlay(ilayer)/1000.E+0
687	ENDDO
688	call profile(nlayer+1,tmp1,tmp2)
689
690	tsurf(1)=tmp2(0)
691	DO ilayer=1,nlayer
692	temp(ilayer)=tmp2(ilayer)
693	ENDDO
694	print*,"check"
695	PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
696	PRINT*,"INPUT TEMPERATURE PROFILE",temp
697
698	c Initialisation albedo / inertie du sol
699	c --------------------------------------
700	albedodat(1)=0.2 ! default value for albedodat
701	PRINT *,'Albedo of bare ground ?'
702	call getin("albedo",albedodat(1))
703	write(,) " albedo = ",albedodat(1)
704
705	inertiedat(1,1)=400 ! default value for inertiedat
706	PRINT *,'Soil thermal inertia (SI) ?'
707	call getin("inertia",inertiedat(1,1))
708	write(,) " inertia = ",inertiedat(1,1)
709
710	! Initialize soil properties and temperature
711	! ------------------------------------------
712	volcapa=1.e6 ! volumetric heat capacity
713	DO isoil=1,nsoil
714	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
715	tsoil(isoil)=tsurf(1) ! soil temperature
716	ENDDO
717
718	! Initialize depths
719	! -----------------
720	do isoil=0,nsoil-1
721	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
722	enddo
723	do isoil=1,nsoil
724	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
725	enddo
726
727	! Initialize methane surface tank
728	! -------------------------------
729	tankCH4=2. ! default value for tankCH4
730
731	! Initialize tracers
732	! ------------------
733	if(tracer) then
734	call initracer2(nq,tname) ! We need tracers names for physdem1 and also to read profiles
735
736	! Inquire if we have tracers input profiles
737
738	do iq=1,nq
739
740	file_prof = "profile_"//trim(noms(iq))
741
742	INQUIRE(file=trim(file_prof),exist=findprof)
743
744	IF (findprof) THEN
745	! Read profile for tracer. First line is surface value.
746	PRINT *, "Input file found for tracer ", noms(iq)
747	OPEN(11,file=trim(file_prof),status='old',form='formatted')
748	PRINT *, "---> reading tracer profile from input file ..."
749	READ(11,*) qsurf(1)
750	DO ilayer=1,nlayer
751	READ (11,*) q(ilayer,iq)
752	ENDDO
753	ELSE
754	PRINT *, "No input profile found for ", trim(noms(iq))
755	PRINT *, "---> we set it to zero !!"
756	ENDIF
757
758	enddo ! nq loop
759
760	endif ! end if tracer
761
762	! Initialize upper atmosphere & chemistry
763	! ---------------------------------------
764
765	! Calculate the # of upper chemistry layers
766	call gr_kim_vervack
767
768	write(,)
769	write(,) " With the compiled vertical grid we found :"
770	write(,) " Number of upper chemistry layers =", nlaykim_up
771
772	if (callchim) then
773
774	if (.NOT.allocated(ykim_up)) then
775	allocate(ykim_up(nkim,1,nlaykim_up))
776	endif
777
778	! Inquire if we have upper chemistry input profiles
779
780	do iq=1,nkim
781
782	file_prof = "profile_"//trim(cnames(iq))//"_up"
783
784	INQUIRE(file=trim(file_prof),exist=findprof)
785
786	IF (findprof) THEN
787	! Read profile for upper chemistry specie.
788	PRINT *, "Input file found for upper chemistry specie ",
789	& trim(cnames(iq))//"_up"
790	OPEN(11,file=trim(file_prof),status='old',form='formatted')
791	PRINT *, "---> reading upper profile from input file ..."
792	DO ilayer=1,nlaykim_up
793	READ (11,*) ykim_up(iq,1,ilayer)
794	ENDDO
795	ELSE
796	PRINT *, "No input profile found for ",
797	& trim(cnames(iq))//"_up"
798	PRINT *, "---> we set it to zero - that's bold !!"
799	ENDIF
800
801	enddo ! nkim loop
802
803	endif ! of if callchim
804
805
806	c Write a "startfi" file
807	c --------------------
808	c This file will be read during the first call to "physiq".
809	c It is needed to transfert physics variables to "physiq"...
810
811	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
812	& dtphys,real(day0),time,cell_area,
813	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
814	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
815	& dtphys,time,
816	& tsurf,tsoil,emis,q2,qsurf,tankCH4)
817
818	c=======================================================================
819	c BOUCLE TEMPORELLE DU MODELE 1D
820	c=======================================================================
821
822	firstcall=.true.
823	lastcall=.false.
824
825	DO idt=1,ndt
826	IF (idt.eq.ndt) then !test
827	lastcall=.true.
828	call stellarlong(day*1.0,zls)
829	! write(103,) 'Ls=',zls180./pi
830	! write(103,) 'Lat=', latitude(1)180./pi
831	! write(103,*) 'RunEnd - Atmos. Temp. File'
832	! write(103,*) 'RunEnd - Atmos. Temp. File'
833	! write(104,) 'Ls=',zls180./pi
834	! write(104,*) 'Lat=', latitude(1)
835	! write(104,*) 'RunEnd - Atmos. Temp. File'
836	ENDIF
837
838	c calcul du geopotentiel
839	c ~~~~~~~~~~~~~~~~~~~~~
840
841
842	DO ilayer=1,nlayer
843
844	! if(autozlevs)then
845	! s(ilayer)=(play(ilayer)/psurf)**rcp
846	! else
847	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
848	! endif
849	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
850	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
851	ENDDO
852
853	! DO ilayer=1,nlayer
854	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
855	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
856	! ENDDO
857	phi(1)=pksh(1)(1.E+0-s(1))
858	DO ilayer=2,nlayer
859	phi(ilayer)=phi(ilayer-1)+
860	& pks(h(ilayer-1)+h(ilayer)).5E+0
861	& *(s(ilayer-1)-s(ilayer))
862
863	ENDDO
864
865	c appel de la physique
866	c --------------------
867
868
869	CALL physiq (1,llm,nq,
870	. tname,
871	, firstcall,lastcall,
872	, day,time,dtphys,
873	, plev,play,phi,
874	, u, v,temp, q,
875	, w,
876	C - sorties
877	s du, dv, dtemp, dq,dpsurf)
878
879
880	c evolution du vent : modele 1D
881	c -----------------------------
882
883	c la physique calcule les derivees temporelles de u et v.
884	c on y rajoute betement un effet Coriolis.
885	c
886	c DO ilayer=1,nlayer
887	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
888	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
889	c ENDDO
890
891	c Pour certain test : pas de coriolis a l'equateur
892	c if(latitude(1).eq.0.) then
893	DO ilayer=1,nlayer
894	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
895	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
896	ENDDO
897	c end if
898	c
899	c
900	c Calcul du temps au pas de temps suivant
901	c ---------------------------------------
902	firstcall=.false.
903	time=time+dtphys/daysec
904	IF (time.gt.1.E+0) then
905	time=time-1.E+0
906	day=day+1
907	ENDIF
908
909	c calcul des vitesses et temperature au pas de temps suivant
910	c ----------------------------------------------------------
911
912	DO ilayer=1,nlayer
913	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
914	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
915	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
916	ENDDO
917
918	c calcul des pressions au pas de temps suivant
919	c ----------------------------------------------------------
920
921	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
922	DO ilevel=1,nlevel
923	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
924	ENDDO
925	DO ilayer=1,nlayer
926	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
927	ENDDO
928
929	c calcul traceur au pas de temps suivant
930	c --------------------------------------
931
932	DO iq = 1, nq
933	DO ilayer=1,nlayer
934	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
935	ENDDO
936	END DO
937
938	c ========================================================
939	c GESTION DES SORTIE
940	c ========================================================
941	if(saveprofile)then
942	OPEN(12,file='profile.out',form='formatted')
943	write(12,*) tsurf
944	DO ilayer=1,llm
945	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
946	ENDDO
947	CLOSE(12)
948	endif
949
950
951	ENDDO ! fin de la boucle temporelle
952
953	write(,) "rcm1d: Everything is cool."
954
955	c ========================================================
956	end !rcm1d
957
958

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

Download in other formats:

Original Format