Context Navigation

rcm1d.F @ 3374

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

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

Download in other formats:

Original Format