Context Navigation

rcm1d.F @ 3955

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

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

Download in other formats:

Original Format