Context Navigation

rcm1d.F @ 3691

Last change on this file since 3691 was 3672, checked in by afalco, 11 months ago
generic/pluto: tsurf/tsoil taken from phys_state_var_mod rather than declared locally so they are updated in the outputs (profile.out, restart1D.nc). AF
Property svn:executable set to ``*
File size: 35.4 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
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, 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, 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	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
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	iphysiq=1 ! in 1D model physics are called evry time step
531	ecritphy=day_step ! default value for ecritphy
532	PRINT *,'Nunber of steps between writediagfi ?'
533	call getin("ecritphy",ecritphy)
534	write(,) " ecritphy = ",ecritphy
535
536	ndt=10 ! default value for ndt
537	PRINT *,'Number of sols to run ?'
538	call getin("ndt",ndt)
539	write(,) " ndt = ",ndt
540	nday=ndt
541
542	ndt=ndt*day_step
543	dtphys=daysec/day_step
544	write(,)"-------------------------------------"
545	write(,)"-------------------------------------"
546	write(,)"--> Physical timestep is ",dtphys
547	write(,)"-------------------------------------"
548	write(,)"-------------------------------------"
549
550	! initializations, as with iniphysiq.F90 for the 3D GCM
551	call init_physics_distribution(regular_lonlat,4,
552	& 1,1,1,nlayer,1)
553	call init_interface_dyn_phys
554	CALL init_regular_lonlat(1,1,longitude,latitude,
555	& (/0.,0./),(/0.,0./))
556	call init_geometry(1,longitude,latitude,
557	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
558	& cell_area)
559	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
560	! call init_vertical_layers(nlayer,preff,scaleheight,
561	! & ap,bp,aps,bps,presnivs,pseudoalt)
562	! call init_dimphy(1,nlayer) ! Initialize dimphy module
563	call ini_planete_mod(nlayer,preff,ap,bp)
564
565	!!! CALL INIFIS
566	!!! - read callphys.def
567	!!! - calculate sine and cosine of longitude and latitude
568	!!! - calculate mugaz and cp
569	!!! NB: some operations are being done dummily in inifis in 1D
570	!!! - physical constants: nevermind, things are done allright below
571	!!! - physical frequency: nevermind, in inifis this is a simple print
572	cpp=-9999. ! dummy init for inifis, will be rewrite later on
573	r=-9999. ! dummy init for inifis, will be rewrite later on
574	CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
575	. latitude,longitude,cell_area,rad,g,r,cpp)
576
577	nsoil=nsoilmx
578
579	! At this point, both getin() and getin_p() functions have been used,
580	! and the run.def file can be removed.
581	call system("rm -rf run.def")
582
583	!!! We check everything is OK.
584	PRINT *,"CHECK"
585	PRINT *,"--> mugaz = ",mugaz
586	PRINT *,"--> cpp = ",cpp
587	r = 8.314511E+0 * 1000.E+0 / mugaz
588	rcp = r / cpp
589
590	c output spectrum?
591	write(,)"Output spectral OLR?"
592	specOLR=.false.
593	call getin("specOLR",specOLR)
594	write(,)" specOLR = ",specOLR
595
596	c option for temperature evolution?
597	write(,)"accelerate temperature?"
598	accelerate_temperature=.false.
599	call getin("accelerate_temperature",accelerate_temperature)
600	write(,)" accelerate_temperature = ",accelerate_temperature
601
602	write(,)"factor for temperature acceleration"
603	factor_acceleration=1.0
604	call getin("factor_acceleration",factor_acceleration)
605	write(,)" factor_acceleration = ",factor_acceleration
606
607
608
609	c
610	c pour le schema d'ondes de gravite
611	c ---------------------------------
612	zmea(1)=0.E+0
613	zstd(1)=0.E+0
614	zsig(1)=0.E+0
615	zgam(1)=0.E+0
616	zthe(1)=0.E+0
617
618	c Initialisation des traceurs
619	c ---------------------------
620
621	DO iq=1,nq
622	DO ilayer=1,nlayer
623	q(ilayer,iq) = 0.
624	ENDDO
625	ENDDO
626
627	DO iq=1,nq
628	qsurf(iq) = 0.
629	ENDDO
630
631	if (tracer) then ! Initialize tracers here.
632
633	write(,) "rcm1d : initializing tracers profiles"
634
635	do iq=1,nq
636
637	if (iq.eq.i_n2) then
638	DO ilayer=1,nlayer
639	q(ilayer,iq) = 1
640	ENDDO
641	else if (iq.eq.i_ch4_gas) then
642	ch4ref=0.5 ! default value for vmr ch4
643	PRINT *,'vmr CH4 (%) ?'
644	call getin("ch4ref",ch4ref)
645	write(,) " CH4 ref (%) = ",ch4ref
646	DO ilayer=1,nlayer
647	q(ilayer,iq) = 0.01ch4ref16./28.
648	ENDDO
649	else if (iq.eq.i_co_gas) then
650	coref=0.05 ! default value for vmr ch4
651	PRINT *,'vmr CO (%) ?'
652	call getin("coref",coref)
653	write(,) " CO ref (%) = ",coref
654	DO ilayer=1,nlayer
655	q(ilayer,iq) = 0.01coref28./28.
656	ENDDO
657	else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30)
658	& .or.(iq.eq.i_haze).or.(iq.eq.i_haze50)
659	& .or.(iq.eq.i_haze100)) then
660	hazeref=0. ! default value for haze mmr
661	PRINT *,'qhaze (kg/kg) ?'
662	call getin("hazeref",hazeref)
663	write(,) " haze ref (kg/kg) = ",hazeref
664	DO ilayer=1,nlayer
665	q(ilayer,iq) = hazeref
666	ENDDO
667	else if (iq.eq.i_prec_haze) then
668	prechazeref=0. ! default value for vmr ch4
669	PRINT *,'qprechaze (kg/kg) ?'
670	call getin("prechazeref",prechazeref)
671	write(,) " prec haze ref (kg/kg) = ",prechazeref
672	DO ilayer=1,nlayer
673	q(ilayer,iq) = prechazeref
674	ENDDO
675	else
676	DO ilayer=1,nlayer
677	q(ilayer,iq) = 0.
678	ENDDO
679	endif
680	enddo ! of do iq=1,nq
681
682	endif ! of tracer
683
684	c Initialisation pour prendre en compte les vents en 1-D
685	c ------------------------------------------------------
686	ptif=2.E+0omegsinlat(1)
687
688
689	c vent geostrophique
690	gru=10. ! default value for gru
691	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
692	call getin("u",gru)
693	write(,) " u = ",gru
694	grv=0. !default value for grv
695	PRINT *,'meridional northward component of the geostrophic',
696	&' wind (m/s) ?'
697	call getin("v",grv)
698	write(,) " v = ",grv
699
700	! To be clean, also set vertical winds to zero
701	w(1:nlayer)=0
702
703	c Initialisation des vents au premier pas de temps
704	DO ilayer=1,nlayer
705	u(ilayer)=gru
706	v(ilayer)=grv
707	ENDDO
708
709	c energie cinetique turbulente
710	DO ilevel=1,nlevel
711	q2(ilevel)=0.E+0
712	ENDDO
713
714	c emissivity / surface n2 ice ( + h2o ice??)
715	c -------------------------------------------
716	emis(1)=emissiv ! default value for emissivity
717	PRINT *,'Emissivity of bare ground ?'
718	call getin("emis",emis(1))
719	write(,) " emis = ",emis(1)
720	emissiv=emis(1) ! we do this so that condense_n2 sets things to the right
721	! value if there is no snow
722
723	if(i_n2.gt.0)then
724	qsurf(i_n2)=0 ! default value for n2ice
725	print*,'Initial n2 ice on the surface (kg.m-2)'
726	call getin("n2ice",qsurf(i_n2))
727	write(,) " n2ice = ",qsurf(i_n2)
728	IF (qsurf(i_n2).ge.1.E+0) THEN
729	! if we have some n2 ice on the surface, change emissivity
730	if (latitude(1).ge.0) then ! northern hemisphere
731	emis(1)=emisice(1)
732	else ! southern hemisphere
733	emis(1)=emisice(2)
734	endif
735	ENDIF
736	endif
737
738	c calcul des pressions et altitudes en utilisant les niveaux sigma
739	c ----------------------------------------------------------------
740
741	c Vertical Coordinates
742	c """"""""""""""""""""
743	hybrid=.true.
744	PRINT *,'Hybrid coordinates ?'
745	call getin("hybrid",hybrid)
746	write(,) " hybrid = ", hybrid
747
748
749	autozlevs=.false.
750	PRINT *,'Auto-discretise vertical levels ?'
751	call getin("autozlevs",autozlevs)
752	write(,) " autozlevs = ", autozlevs
753
754	pceil = psurf / 1000.0 ! Pascals
755	PRINT *,'Ceiling pressure (Pa) ?'
756	call getin("pceil",pceil)
757	write(,) " pceil = ", pceil
758
759	! Test of incompatibility:
760	! if autozlevs used, cannot have hybrid too
761
762	if (autozlevs.and.hybrid) then
763	print*,'Cannot use autozlevs and hybrid together!'
764	call abort
765	endif
766
767	if(autozlevs)then
768
769	open(91,file="z2sig.def",form='formatted')
770	read(91,*) Hscale
771	DO ilayer=1,nlayer-2
772	read(91,*) Hmax
773	enddo
774	close(91)
775
776
777	print*,'Hmax = ',Hmax,' km'
778	print*,'Auto-shifting Hscale to:'
779	! Hscale = Hmax / log(psurf/100.0)
780	Hscale = Hmax / log(psurf/pceil)
781	print*,'Hscale = ',Hscale,' km'
782
783	! none of this matters if we dont care about zlay
784
785	endif
786
787	call disvert_noterre
788	! now that disvert has been called, initialize module vertical_layers_mod
789	call init_vertical_layers(nlayer,preff,scaleheight,
790	& ap,bp,aps,bps,presnivs,pseudoalt)
791
792	if(.not.autozlevs)then
793	! we want only the scale height from z2sig, in order to compute zlay
794	open(91,file="z2sig.def",form='formatted')
795	read(91,*) Hscale
796	close(91)
797	endif
798
799	! if(autozlevs)then
800	! open(94,file="Hscale.temp",form='formatted')
801	! read(94,*) Hscale
802	! close(94)
803	! endif
804
805	DO ilevel=1,nlevel
806	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
807	ENDDO
808
809	DO ilayer=1,nlayer
810	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
811	ENDDO
812
813
814
815	DO ilayer=1,nlayer
816	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
817	! & /g
818	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
819	ENDDO
820
821	! endif
822
823	c profil de temperature au premier appel
824	c --------------------------------------
825	pks=psurf**rcp
826
827	c altitude en km dans profile: on divise zlay par 1000
828	tmp1(0)=0.E+0
829	DO ilayer=1,nlayer
830	tmp1(ilayer)=zlay(ilayer)/1000.E+0
831	ENDDO
832	call profile(nlayer+1,tmp1,tmp2)
833
834	tsurf(1)=tmp2(0)
835	DO ilayer=1,nlayer
836	temp(ilayer)=tmp2(ilayer)
837	ENDDO
838	print*,"check"
839	PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
840	PRINT*,"INPUT TEMPERATURE PROFILE",temp
841
842	c Initialisation albedo / inertie du sol
843	c --------------------------------------
844	albedodat(1)=0.2 ! default value for albedodat
845	PRINT *,'Albedo of bare ground ?'
846	call getin("albedo",albedodat(1))
847	write(,) " albedo = ",albedodat(1)
848
849	inertiedat(1,1)=400 ! default value for inertiedat
850	PRINT *,'Soil thermal inertia (SI) ?'
851	call getin("inertia",inertiedat(1,1))
852	write(,) " inertia = ",inertiedat(1,1)
853
854	! Initialize soil properties and temperature
855	! ------------------------------------------
856	volcapa=1.e6 ! volumetric heat capacity
857	lecttsoil=0 ! default value for lecttsoil
858	call getin("lecttsoil",lecttsoil)
859	if (lecttsoil==1) then
860	OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
861	DO isoil=1,nsoil
862	READ (14,*) tsoil(1,isoil)
863	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
864	ENDDO
865	GOTO 201
866	101 STOP'fichier proftsoil inexistant'
867	201 CONTINUE
868	CLOSE(14)
869
870	else
871	DO isoil=1,nsoil
872	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
873	tsoil(1,isoil)=tsurf(1) ! soil temperature
874	ENDDO
875	endif
876
877
878	! Initialize depths
879	! -----------------
880	do isoil=0,nsoil-1
881	mlayer(isoil)=3.e-2(2.*(isoil-0.5)) ! mid-layer depth
882	enddo
883	do isoil=1,nsoil
884	layer(isoil)=3.e-2(2.*(isoil-1)) ! layer depth
885	enddo
886
887	! Initialize haze profile
888	! ------------------------------------------
889	if (haze) then
890
891	lecthaze=0 ! default value for lecthaze
892	call getin("lecthaze",lecthaze)
893	if (lecthaze==1) then
894	OPEN(15,file='profhaze',status='old',form='formatted',err=301)
895	DO iq=1,nq
896	if (iq.eq.i_haze) then
897	DO ilayer=1,nlayer
898	READ (15,*) q(ilayer,iq)
899	ENDDO
900	endif
901	ENDDO
902	GOTO 401
903	301 STOP'Problem with profhaze file'
904	401 CONTINUE
905	CLOSE(15)
906	endif
907	endif
908	! Initialize cloud fraction and oceanic ice !AF24: removed
909	! Initialize slab ocean !AF24: removed
910	! Initialize chemical species !AF24: removed photochem
911
912
913	c Write a "startfi" file
914	c --------------------
915	c This file will be read during the first call to "physiq".
916	c It is needed to transfert physics variables to "physiq"...
917
918	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
919	& dtphys,real(day0),time,cell_area,
920	& albedodat,zmea,zstd,zsig,zgam,zthe)
921	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
922	& dtphys,time,
923	& tsurf,tsoil,inertiedat,emis,albedodat,q2,qsurf)
924
925	c=======================================================================
926	c BOUCLE TEMPORELLE DU MODELE 1D
927	c=======================================================================
928
929	firstcall=.true.
930	lastcall=.false.
931
932	DO idt=1,ndt
933	IF (idt.eq.ndt) then !test
934	lastcall=.true.
935	call stellarlong(day*1.0,zls)
936	! write(103,) 'Ls=',zls180./pi
937	! write(103,) 'Lat=', latitude(1)180./pi
938	! write(103,*) 'RunEnd - Atmos. Temp. File'
939	! write(103,*) 'RunEnd - Atmos. Temp. File'
940	! write(104,) 'Ls=',zls180./pi
941	! write(104,*) 'Lat=', latitude(1)
942	! write(104,*) 'RunEnd - Atmos. Temp. File'
943	ENDIF
944
945	c calcul du geopotentiel
946	c ~~~~~~~~~~~~~~~~~~~~~
947
948
949	DO ilayer=1,nlayer
950
951	! if(autozlevs)then
952	! s(ilayer)=(play(ilayer)/psurf)**rcp
953	! else
954	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
955	! endif
956	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
957	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
958	ENDDO
959
960	! DO ilayer=1,nlayer
961	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
962	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
963	! ENDDO
964	phi(1)=pksh(1)(1.E+0-s(1))
965	DO ilayer=2,nlayer
966	phi(ilayer)=phi(ilayer-1)+
967	& pks(h(ilayer-1)+h(ilayer)).5E+0
968	& *(s(ilayer-1)-s(ilayer))
969
970	ENDDO
971
972	c appel de la physique
973	c --------------------
974
975
976	CALL physiq (1,llm,nq,
977	, firstcall,lastcall,
978	, day,time,dtphys,
979	, plev,play,phi,
980	, u, v,temp, q,
981	, w,
982	C - sorties
983	s du, dv, dtemp, dq,dpsurf)
984
985
986	c evolution du vent : modele 1D
987	c -----------------------------
988
989	c la physique calcule les derivees temporelles de u et v.
990	c on y rajoute betement un effet Coriolis.
991	c
992	c DO ilayer=1,nlayer
993	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
994	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
995	c ENDDO
996
997	c Pour certain test : pas de coriolis a l'equateur
998	c if(latitude(1).eq.0.) then
999	DO ilayer=1,nlayer
1000	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
1001	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
1002	ENDDO
1003	c end if
1004	c
1005	c
1006	c Calcul du temps au pas de temps suivant
1007	c ---------------------------------------
1008	firstcall=.false.
1009	time=time+dtphys/daysec
1010	IF (time.gt.1.E+0) then
1011	time=time-1.E+0
1012	day=day+1
1013	ENDIF
1014
1015	c calcul des vitesses et temperature au pas de temps suivant
1016	c ----------------------------------------------------------
1017
1018	DO ilayer=1,nlayer
1019	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
1020	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
1021	if(accelerate_temperature) then
1022	temp(ilayer)=temp(ilayer)+dtphysdtemp(ilayer)
1023	& max(1.0,play(ilayer)1e-5)*factor_acceleration
1024	else
1025	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
1026	endif
1027	ENDDO
1028
1029	c calcul des pressions au pas de temps suivant
1030	c ----------------------------------------------------------
1031
1032	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
1033	DO ilevel=1,nlevel
1034	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
1035	ENDDO
1036	DO ilayer=1,nlayer
1037	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
1038	ENDDO
1039
1040	c calcul traceur au pas de temps suivant
1041	c --------------------------------------
1042
1043	DO iq = 1, nq
1044	DO ilayer=1,nlayer
1045	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
1046	ENDDO
1047	END DO
1048
1049	c ========================================================
1050	c GESTION DES SORTIE
1051	c ========================================================
1052	if(saveprofile)then
1053	OPEN(12,file='profile.out',form='formatted')
1054	write(12,*) tsurf
1055	DO ilayer=1,llm
1056	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
1057	ENDDO
1058	CLOSE(12)
1059	endif
1060	! save haze profile
1061	if (haze.and.lecthaze.eq.1) then
1062	OPEN(16,file='profhaze.out',form='formatted')
1063	DO iq=1,nq
1064	if (iq.eq.i_haze) then
1065	DO ilayer=1,nlayer
1066	write(16,*) q(ilayer,iq)
1067	ENDDO
1068	endif
1069	ENDDO
1070	CLOSE(16)
1071	endif
1072
1073	ENDDO ! fin de la boucle temporelle
1074
1075	write(,) "rcm1d: Everything is cool."
1076
1077	c ========================================================
1078	end !rcm1d
1079
1080

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

Download in other formats:

Original Format