Context Navigation

rcm1d.F @ 2743

Last change on this file since 2743 was 2366, checked in by jvatant, 6 years ago

Titan GCM : Major maintenance catching up commits from the generic including :

r2356 and 2354 removing obsolete old dynamical core
various minor addition to physics and gestion of phys_state_var_mode, especially in dyn1d
adding MESOSCALE CPP keys around chemistry and microphysics (disabled in mesoscale for now)

Property svn:executable set to *

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

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

Download in other formats:

Original Format