Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

rcm1d.F

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

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

Download in other formats:

Original Format