Context Navigation

rcm1d.F @ 3772

Last change on this file since 3772 was 3758, checked in by afalco, 11 months ago
Pluto: rcm1d: do not import iphysiq (not existing anymore) Make Lott scheme deactivated by default. 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, 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 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 during the run
83	INTEGER day_step ! number of time steps per day
84	REAL time ! time (0<time<1 ; time=0.5 a midi)
85	REAL play(llm) ! Pressure at the middle of the layers (Pa)
86	REAL plev(llm+1) ! intermediate pressure levels (pa)
87	REAL psurf
88	REAL u(llm),v(llm) ! zonal, meridional wind
89	REAL gru,grv ! prescribed "geostrophic" background wind
90	REAL temp(llm) ! temperature at the middle of the layers
91	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
92	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
93	! REAL n2ice ! n2ice layer (kg.m-2) !not used anymore
94	integer :: i_n2=0 ! tracer index of n2 ice
95	integer :: i_ch4_ice=0 ! tracer index of ch4 ice
96	integer :: i_ch4_gas=0 ! tracer index of ch4 gas
97	integer :: i_co_ice=0 ! tracer index of co ice
98	integer :: i_co_gas=0 ! tracer index of co gas
99	integer :: i_prec_haze=0 ! tracer index of haze
100	integer :: i_haze=0 ! tracer index of haze
101	integer :: i_haze10=0 ! tracer index of haze
102	integer :: i_haze30=0 ! tracer index of haze
103	integer :: i_haze50=0 ! tracer index of haze
104	integer :: i_haze100=0 ! tracer index of haze
105
106	REAL emis(1) ! surface layer
107	REAL q2(llm+1) ! Turbulent Kinetic Energy
108	REAL zlay(llm) ! altitude estimee dans les couches (km)
109
110	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
111	REAL du(llm),dv(llm),dtemp(llm)
112	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
113	REAL dpsurf(1)
114	REAL,ALLOCATABLE :: dq(:,:)
115	REAL,ALLOCATABLE :: dqdyn(:,:)
116
117	c Various intermediate variables
118	! INTEGER thermo
119	REAL zls
120	REAL phi(llm),h(llm),s(llm)
121	REAL pks, ptif, w(llm)
122	INTEGER ierr, aslun
123	REAL tmp1(0:llm),tmp2(0:llm)
124	integer :: nq !=1 ! number of tracers
125
126	character*2 str2
127	character (len=7) :: str7
128	character(len=44) :: txt
129	character(len=44) :: tracer_profile_file_name
130
131	logical oldcompare, earthhack,saveprofile
132
133	! added by RW for zlay computation
134	real Hscale, Hmax, rho, dz
135
136	! pluto specific
137	real ch4ref ! % ch4
138	real coref ! % co
139	real hazeref ! haze kg/kg
140	real prechazeref ! prec haze kg/kg
141
142
143	! added by RW for autozlevs computation
144	logical autozlevs
145	real nu, xx, pMIN, zlev, Htop
146	real logplevs(llm)
147
148	! added by BC to accelerate convergence
149	logical accelerate_temperature
150	real factor_acceleration
151
152	! added by AS to avoid the use of adv trac common
153	character*30,allocatable :: nametmp(:) !
154
155	real :: latitude(1), longitude(1), cell_area(1)
156
157	! added by JVO and YJ to read modern traceur.def
158	character(len=500) :: line ! to store a line of text
159	LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def
160
161	c=======================================================================
162	c INITIALISATION
163	c=======================================================================
164	! check if 'rcm1d.def' file is around
165	open(90,file='rcm1d.def',status='old',form='formatted',
166	& iostat=ierr)
167	if (ierr.ne.0) then
168	write(,) 'Cannot find required file "rcm1d.def"'
169	write(,) 'which should contain some input parameters'
170	write(,) ' ... might as well stop here ...'
171	stop
172	else
173	close(90)
174	endif
175
176	! now, run.def is needed anyway. so we create a dummy temporary one
177	! for ioipsl to work. if a run.def is already here, stop the
178	! program and ask the user to do a bit of cleaning
179	open(90,file='run.def',status='old',form='formatted',
180	& iostat=ierr)
181	if (ierr.eq.0) then
182	close(90)
183	write(,) 'There is already a run.def file.'
184	write(,) 'This is not compatible with 1D runs.'
185	write(,) 'Please remove the file and restart the run.'
186	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
187	stop
188	else
189	call system('touch run.def')
190	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
191	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
192	endif
193
194	! read nq from traceur.def
195	open(90,file='traceur.def',status='old',form='formatted',
196	& iostat=ierr)
197	if (ierr.eq.0) then
198	! read number of tracers:
199	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
200	READ(90,'(A)') line
201	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
202	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
203	ELSE
204	moderntracdef = .true.
205	DO
206	READ(90,'(A)',iostat=ierr) line
207	IF (ierr.eq.0) THEN
208	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
209	READ(line,*,iostat=ierr) nq
210	EXIT
211	ENDIF
212	ELSE ! If pb, or if reached EOF without having found nbtr
213	write(,) "rcm1d: error reading number of tracers"
214	write(,) " (first line of traceur.def) "
215	stop
216	ENDIF
217	ENDDO
218	ENDIF ! if modern or standard traceur.def
219	else
220	nq=0
221	endif
222	close(90)
223
224	! Initialize dimphy module
225	call init_dimphy(1,llm)
226
227	! now initialize arrays using phys_state_var_init
228	! but first initialise naerkind (from callphys.def)
229	naerkind=0 !default
230	call getin("naerkind",naerkind)
231
232	call phys_state_var_init(nq)
233
234	saveprofile=.false.
235	saveprofile=.true.
236
237	c ----------------------------------------
238	c Default values (corresponding to Mars)
239	c ----------------------------------------
240
241	pi=2.E+0*asin(1.E+0)
242
243	c Parametres Couche limite et Turbulence
244	c --------------------------------------
245	z0 = 1.e-2 ! surface roughness (m) ~0.01
246	emin_turb = 1.e-6 ! energie minimale ~1.e-8
247	lmixmin = 30 ! longueur de melange ~100
248
249	c propriete optiques des calottes et emissivite du sol
250	c ----------------------------------------------------
251	emissiv= 0.95 ! Emissivite du sol martien ~.95
252	emisice(1)=0.95 ! Emissivite calotte nord
253	emisice(2)=0.95 ! Emissivite calotte sud
254
255	iceradius(1) = 100.e-6 ! mean scat radius of n2 snow (north)
256	iceradius(2) = 100.e-6 ! mean scat radius of n2 snow (south)
257	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
258	dtemisice(2) = 2. ! time scale for snow metamorphism (south
259	hybrid=.false.
260
261	c ------------------------------------------------------
262	c Load parameters from "run.def" and "gases.def"
263	c ------------------------------------------------------
264
265
266	! check if we are going to run with or without tracers
267	write(,) "Run with or without tracer transport ?"
268	tracer=.false. ! default value
269	call getin("tracer",tracer)
270	write(,) " tracer = ",tracer
271
272	! OK. now that run.def has been read once -- any variable is in memory.
273	! so we can dump the dummy run.def
274	! call system("rm -rf run.def") ! Ehouarn: delay this to after inifis
275
276	! while we're at it, check if there is a 'traceur.def' file
277	! and preocess it, if necessary. Otherwise initialize tracer names
278	if (tracer) then
279	! load tracer names from file 'traceur.def'
280	open(90,file='traceur.def',status='old',form='formatted',
281	& iostat=ierr)
282	if (ierr.eq.0) then
283	write(,) "rcm1d: Reading file traceur.def"
284	! read number of tracers:
285	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
286	READ(90,'(A)') line
287	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
288	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
289	ELSE
290	moderntracdef = .true.
291	DO
292	READ(90,'(A)',iostat=ierr) line
293	IF (ierr.eq.0) THEN
294	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
295	READ(line,*,iostat=ierr) nq
296	EXIT
297	ENDIF
298	ELSE ! If pb, or if reached EOF without having found nbtr
299	write(,) "rcm1d: error reading number of tracers"
300	write(,) " (first line of traceur.def) "
301	stop
302	ENDIF
303	ENDDO
304	ENDIF ! if modern or standard traceur.def
305	nqtot=nq ! set value of nqtot (in infotrac module) as nq
306	if (ierr.ne.0) then
307	write(,) "rcm1d: error reading number of tracers"
308	write(,) " (first line of traceur.def) "
309	stop
310	endif
311	if (nq>0) then
312	allocate(tname(nq))
313	allocate(noms(nq))
314	allocate(q(llm,nq))
315	allocate(qsurf(nq))
316	allocate(dq(llm,nq))
317	allocate(dqdyn(llm,nq))
318	else
319	write(,) "rcm1d: Error. You set tracer=.true."
320	write(,) " but # of tracers in traceur.def is ",nq
321	stop
322	endif
323
324	do iq=1,nq
325	! minimal version, just read in the tracer names, 1 per line
326	read(90,*,iostat=ierr) tname(iq)
327	noms(iq)=tname(iq)
328	if (ierr.ne.0) then
329	write(,) 'rcm1d: error reading tracer names...'
330	stop
331	endif
332	enddo !of do iq=1,nq
333	! check for n2 / ch4_ice tracers:
334	i_n2=0
335	do iq=1,nq
336	if (tname(iq)=="n2") then
337	i_n2=iq
338	elseif (tname(iq)=="ch4_ice") then
339	i_ch4_ice=iq
340	elseif (tname(iq)=="co_ice") then
341	i_co_ice=iq
342	elseif (tname(iq)=="ch4_gas") then
343	i_ch4_gas=iq
344	elseif (tname(iq)=="co_gas") then
345	i_co_gas=iq
346	elseif (tname(iq)=="haze") then
347	i_haze=iq
348	endif
349	enddo
350	else
351	write(,) 'Cannot find required file "traceur.def"'
352	write(,) ' If you want to run with tracers, I need it'
353	write(,) ' ... might as well stop here ...'
354	stop
355	endif
356	close(90)
357
358
359	else ! of if (tracer)
360	nqtot=0
361	nq=0
362	! still, make allocations for 1 dummy tracer
363	allocate(tname(1))
364	allocate(qsurf(1))
365	allocate(q(llm,1))
366	allocate(dq(llm,1))
367
368	! Check that tracer boolean is compliant with number of tracers
369	! -- otherwise there is an error (and more generally we have to be consistent)
370	if (nq .ge. 1) then
371	write(,) "------------------------------"
372	write(,) "rcm1d: You set tracer=.false."
373	write(,) " But set number of tracers to ",nq
374	write(,) " > If you want tracers, set tracer=.true."
375	write(,) "------------------------------"
376	stop
377	endif
378	endif ! of if (tracer)
379
380	!!! GEOGRAPHICAL INITIALIZATIONS
381	!!! AREA. useless in 1D
382	cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
383	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
384	phisfi(1)=0.E+0
385	!!! LATITUDE. can be set.
386	latitude=0 ! default value for latitude
387	PRINT *,'latitude (in degrees) ?'
388	call getin("latitude",latitude)
389	write(,) " latitude = ",latitude
390	latitude=latitude*pi/180.E+0
391	!!! LONGITUDE. useless in 1D.
392	longitude=0.E+0
393	longitude=longitude*pi/180.E+0
394
395
396	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397	!!!! PLANETARY CONSTANTS !!!!
398	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399
400	g = -99999.
401	PRINT *,'GRAVITY in m s-2 ?'
402	call getin("g",g)
403	IF (g.eq.-99999.) THEN
404	PRINT *,"STOP. I NEED g IN RCM1D.DEF."
405	STOP
406	ELSE
407	PRINT *,"--> g = ",g
408	ENDIF
409
410	rad = -99999.
411	PRINT *,'PLANETARY RADIUS in m ?'
412	call getin("rad",rad)
413	! Planetary radius is needed to compute shadow of the rings
414	IF (rad.eq.-99999.) THEN
415	PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
416	STOP
417	ELSE
418	PRINT *,"--> rad = ",rad
419	ENDIF
420
421	daysec = -99999.
422	PRINT *,'LENGTH OF A DAY in s ?'
423	call getin("daysec",daysec)
424	IF (daysec.eq.-99999.) THEN
425	PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
426	STOP
427	ELSE
428	PRINT *,"--> daysec = ",daysec
429	ENDIF
430	omeg=4.*asin(1.)/(daysec)
431	PRINT *,"OK. FROM THIS I WORKED OUT:"
432	PRINT *,"--> omeg = ",omeg
433
434	year_day = -99999.
435	PRINT *,'LENGTH OF A YEAR in days ?'
436	call getin("year_day",year_day)
437	IF (year_day.eq.-99999.) THEN
438	PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
439	STOP
440	ELSE
441	PRINT *,"--> year_day = ",year_day
442	ENDIF
443
444	periastr = -99999.
445	PRINT *,'MIN DIST STAR-PLANET in AU ?'
446	call getin("periastr",periastr)
447	IF (periastr.eq.-99999.) THEN
448	PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
449	STOP
450	ELSE
451	PRINT *,"--> periastr = ",periastr
452	ENDIF
453
454	apoastr = -99999.
455	PRINT *,'MAX DIST STAR-PLANET in AU ?'
456	call getin("apoastr",apoastr)
457	IF (apoastr.eq.-99999.) THEN
458	PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
459	STOP
460	ELSE
461	PRINT *,"--> apoastr = ",apoastr
462	ENDIF
463
464	peri_day = -99999.
465	PRINT *,'DATE OF PERIASTRON in days ?'
466	call getin("peri_day",peri_day)
467	IF (peri_day.eq.-99999.) THEN
468	PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
469	STOP
470	ELSE IF (peri_day.gt.year_day) THEN
471	PRINT *,"STOP. peri_day.gt.year_day"
472	STOP
473	ELSE
474	PRINT *,"--> peri_day = ", peri_day
475	ENDIF
476
477	obliquit = -99999.
478	PRINT *,'OBLIQUITY in deg ?'
479	call getin("obliquit",obliquit)
480	IF (obliquit.eq.-99999.) THEN
481	PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
482	STOP
483	ELSE
484	PRINT *,"--> obliquit = ",obliquit
485	ENDIF
486
487	psurf = -99999.
488	PRINT *,'SURFACE PRESSURE in Pa ?'
489	call getin("psurf",psurf)
490	IF (psurf.eq.-99999.) THEN
491	PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
492	STOP
493	ELSE
494	PRINT *,"--> psurf = ",psurf
495	ENDIF
496	!! we need reference pressures
497	pa=psurf/30.
498	preff=psurf ! these values are not needed in 1D (are you sure JL12)
499
500	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
501	!!!! END PLANETARY CONSTANTS !!!!
502	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503
504	c Date et heure locale du debut du run
505	c ------------------------------------
506	c Date (en sols depuis le solstice de printemps) du debut du run
507	day0 = 0 ! default value for day0
508	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
509	call getin("day0",day0)
510	day=float(day0)
511	write(,) " day0 = ",day0
512	c Heure de demarrage
513	time=0 ! default value for time
514	write(,)'Initial local time (in hours, between 0 and 24)?'
515	call getin("time",time)
516	write(,)" time = ",time
517	time=time/24.E+0 ! convert time (hours) to fraction of sol
518
519
520	c Discretisation (Definition de la grille et des pas de temps)
521	c --------------
522	c
523	nlayer=llm
524	nlevel=nlayer+1
525
526	day_step=48 ! default value for day_step
527	PRINT *,'Number of time steps per sol ?'
528	call getin("day_step",day_step)
529	write(,) " day_step = ",day_step
530
531	diagfi_output_rate=24 ! default value for diagfi_output_rate
532	PRINT *,'Nunber of steps between writediagfi ?'
533	call getin("diagfi_output_rate",diagfi_output_rate)
534	write(,) " diagfi_output_rate = ",diagfi_output_rate
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,[1])
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