Context Navigation

rcm1d.F @ 773

Last change on this file since 773 was 773, checked in by jleconte, 12 years ago

05/09/2012 == JL

Correction of the calculation of the solar longitude in tlocked case.

-Can now handle any prograde resonance with nres=omega_rot/omega_orb.
-Sun now goes westward for the standard 2:1 case, as expected.

In the gray case, the separation between kappa_IR and VI is now set by

wave number, independently of the usual IR/VISIBLE calculation separation.
i.e. kappa_IR can be used in the calculation of the downward stellar flux

if the wavenumber in the band is low enough and vice versa.

In ave_stelspec, stellar flux averaging has been generalized to incorporate

very red/blue stellar spectra (great care must however be taken of the band
limit used for the corralated k distributions).

-Brown dwarf spectra from Allard et al. have been added.
-Any Black body temperature can now be used.

Property svn:executable set to *

File size: 25.9 KB

Line
1	program rcm1d
2
3	! to use 'getin'
4	use ioipsl_getincom
5
6	implicit none
7
8	!==================================================================
9	!
10	! Purpose
11	! -------
12	! Run the physics package of the universal model in a 1D column.
13	!
14	! It can be compiled with a command like (e.g. for 25 layers):
15	! "makegcm -p std -d 25 rcm1d"
16	! It requires the files "callphys.def", "z2sig.def",
17	! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
18	!
19	! Authors
20	! -------
21	! Frederic Hourdin
22	! R. Fournier
23	! F. Forget
24	! F. Montmessin (water ice added)
25	! R. Wordsworth
26	! B. Charnay
27	! A. Spiga
28	! J. Leconte (2012)
29	!
30	!==================================================================
31
32	#include "dimensions.h"
33	#include "paramet.h"
34	#include "dimphys.h"
35	#include "comgeomfi.h"
36	#include "surfdat.h"
37	#include "comsoil.h"
38	#include "comdiurn.h"
39	#include "callkeys.h"
40	#include "comcstfi.h"
41	#include "planete.h"
42	#include "comsaison.h"
43	#include "control.h"
44	#include "comvert.h"
45	#include "netcdf.inc"
46	#include "watercap.h"
47	#include "fisice.h"
48	#include "logic.h"
49	#include "advtrac.h"
50	#include "comgeom.h"
51
52	c --------------------------------------------------------------
53	c Declarations
54	c --------------------------------------------------------------
55	c
56	INTEGER unitstart ! unite d'ecriture de "startfi"
57	INTEGER nlayer,nlevel,nsoil,ndt
58	INTEGER ilayer,ilevel,isoil,idt,iq
59	LOGICAl firstcall,lastcall
60	c
61	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
62	REAL day ! date durant le run
63	REAL time ! time (0<time<1 ; time=0.5 a midi)
64	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
65	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
66	REAL psurf,tsurf
67	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
68	REAL gru,grv ! prescribed "geostrophic" background wind
69	REAL temp(nlayermx) ! temperature at the middle of the layers
70	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
71	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
72	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
73	! REAL co2ice ! co2ice layer (kg.m-2) !not used anymore
74	integer :: i_co2_ice=0 ! tracer index of co2 ice
75	integer :: i_h2o_ice=0 ! tracer index of h2o ice
76	integer :: i_h2o_vap=0 ! tracer index of h2o vapor
77	REAL emis ! surface layer
78	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
79	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
80
81	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
82	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
83	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
84	REAL dpsurf
85	REAL dq(nlayermx,nqmx)
86	REAL dqdyn(nlayermx,nqmx)
87
88	c Various intermediate variables
89	! INTEGER thermo
90	REAL zls
91	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
92	REAL pks, ptif, w(nlayermx)
93	REAL qtotinit, mqtot(nqmx),qtot
94	INTEGER ierr, aslun
95	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
96	Logical tracerdyn
97	integer :: nq=1 ! number of tracers
98
99	character*2 str2
100	character (len=7) :: str7
101
102	logical oldcompare, earthhack,saveprofile
103
104	! added by RW for zlay computation
105	real Hscale, Hmax, rho, dz
106
107	! added by RW for autozlevs computation
108	real nu, xx, pMIN, zlev, Htop
109	real logplevs(nlayermx)
110
111	! added by BC
112	REAL cloudfrac(ngridmx,nlayermx)
113	REAL hice(ngridmx),totcloudfrac(ngridmx)
114	REAL qzero1D !initial water amount on the ground
115
116	c=======================================================================
117	c INITIALISATION
118	c=======================================================================
119
120	saveprofile=.false.
121	saveprofile=.true.
122
123	c ----------------------------------------
124	c Default values (corresponding to Mars)
125	c ----------------------------------------
126
127	pi=2.E+0*asin(1.E+0)
128
129	c Parametres Couche limite et Turbulence
130	c --------------------------------------
131	z0 = 1.e-2 ! surface roughness (m) ~0.01
132	emin_turb = 1.e-6 ! energie minimale ~1.e-8
133	lmixmin = 30 ! longueur de melange ~100
134
135	c propriete optiques des calottes et emissivite du sol
136	c ----------------------------------------------------
137	emissiv= 0.95 ! Emissivite du sol martien ~.95
138	emisice(1)=0.95 ! Emissivite calotte nord
139	emisice(2)=0.95 ! Emissivite calotte sud
140
141	albedice(1)=0.5 ! Albedo calotte nord
142	albedice(2)=0.5 ! Albedo calotte sud
143	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
144	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
145	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
146	dtemisice(2) = 2. ! time scale for snow metamorphism (south
147	hybrid=.false.
148
149	c ------------------------------------------------------
150	c Load parameters from "run.def" and "gases.def"
151	c ------------------------------------------------------
152
153	! check if 'rcm1d.def' file is around
154	open(90,file='rcm1d.def',status='old',form='formatted',
155	& iostat=ierr)
156	if (ierr.ne.0) then
157	write(,) 'Cannot find required file "rcm1d.def"'
158	write(,) 'which should contain some input parameters'
159	write(,) ' ... might as well stop here ...'
160	stop
161	else
162	close(90)
163	endif
164
165	! now, run.def is needed anyway. so we create a dummy temporary one
166	! for ioipsl to work. if a run.def is already here, stop the
167	! program and ask the user to do a bit of cleaning
168	open(90,file='run.def',status='old',form='formatted',
169	& iostat=ierr)
170	if (ierr.eq.0) then
171	close(90)
172	write(,) 'There is already a run.def file.'
173	write(,) 'This is not compatible with 1D runs.'
174	write(,) 'Please remove the file and restart the run.'
175	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
176	stop
177	else
178	call system('touch run.def')
179	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
180	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
181	endif
182
183	! check if we are going to run with or without tracers
184	write(,) "Run with or without tracer transport ?"
185	tracer=.false. ! default value
186	call getin("tracer",tracer)
187	write(,) " tracer = ",tracer
188
189	! OK. now that run.def has been read once -- any variable is in memory.
190	! so we can dump the dummy run.def
191	call system("rm -rf run.def")
192
193	! while we're at it, check if there is a 'traceur.def' file
194	! and preocess it, if necessary. Otherwise initialize tracer names
195	if (tracer) then
196	! load tracer names from file 'traceur.def'
197	open(90,file='traceur.def',status='old',form='formatted',
198	& iostat=ierr)
199	if (ierr.eq.0) then
200	write(,) "rcm1d: Reading file traceur.def"
201	! read number of tracers:
202	read(90,*,iostat=ierr) nq
203	if (ierr.ne.0) then
204	write(,) "rcm1d: error reading number of tracers"
205	write(,) " (first line of traceur.def) "
206	stop
207	else
208	! check that the number of tracers is indeed nqmx
209	if (nq.ne.nqmx) then
210	write(,) "rcm1d: error, wrong number of tracers:"
211	write(,) "nq=",nq," whereas nqmx=",nqmx
212	stop
213	endif
214	endif
215
216	! initialize advection schemes to Van-Leer for all tracers
217	do iq=1,nq
218	iadv(iq)=3 ! Van-Leer
219	enddo
220
221	do iq=1,nq
222	! minimal version, just read in the tracer names, 1 per line
223	read(90,*,iostat=ierr) tnom(iq)
224	if (ierr.ne.0) then
225	write(,) 'rcm1d: error reading tracer names...'
226	stop
227	endif
228	enddo !of do iq=1,nq
229	! check for co2_ice / h2o_ice tracers:
230	i_co2_ice=0
231	i_h2o_ice=0
232	i_h2o_vap=0
233	do iq=1,nq
234	if (tnom(iq)=="co2_ice") then
235	i_co2_ice=iq
236	elseif (tnom(iq)=="h2o_ice") then
237	i_h2o_ice=iq
238	elseif (tnom(iq)=="h2o_vap") then
239	i_h2o_vap=iq
240	endif
241	enddo
242	!if (i_co2_ice==0) then
243	! write(,) "rcm1d: error, we need a 'co2_ice' tracer"
244	! write(,) " (add one to traceur.def)"
245	! stop
246	!endif
247	else
248	write(,) 'Cannot find required file "traceur.def"'
249	write(,) ' If you want to run with tracers, I need it'
250	write(,) ' ... might as well stop here ...'
251	stop
252	endif
253	close(90)
254
255	else
256	! we still need to set (dummy) tracer names for physdem1
257	nq=nqmx
258	do iq=1,nq
259	write(str7,'(a1,i2.2)')'q',iq
260	tnom(iq)=str7
261	enddo
262	! actually, we'll need at least one "co2_ice" tracer
263	! (for surface CO2 ice)
264	tnom(1)="co2_ice"
265	i_co2_ice=1
266	endif ! of if (tracer)
267
268	!!! We have to check that check_cpp_match is F for 1D computations
269	!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
270	check_cpp_match = .false.
271	call getin("check_cpp_match",check_cpp_match)
272	if (check_cpp_match) then
273	print*,"In 1D modeling, check_cpp_match is supposed to be F"
274	print*,"Please correct callphys.def"
275	stop
276	endif
277
278	!!! GEOGRAPHICAL INITIALIZATIONS
279	!!! AREA. useless in 1D
280	area(1)=1.E+0
281	aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
282	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
283	phisfi(1)=0.E+0
284	!!! LATITUDE. can be set.
285	lati(1)=0 ! default value for lati(1)
286	PRINT *,'latitude (in degrees) ?'
287	call getin("latitude",lati(1))
288	write(,) " latitude = ",lati(1)
289	lati(1)=lati(1)*pi/180.E+0
290	!!! LONGITUDE. useless in 1D.
291	long(1)=0.E+0
292	long(1)=long(1)*pi/180.E+0
293
294	!!! CALL INIFIS
295	!!! - read callphys.def
296	!!! - calculate sine and cosine of longitude and latitude
297	!!! - calculate mugaz and cp
298	!!! NB: some operations are being done dummily in inifis in 1D
299	!!! - physical constants: nevermind, things are done allright below
300	!!! - physical frequency: nevermind, in inifis this is a simple print
301	cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
302	CALL inifis(1,llm,day0,daysec,dtphys,
303	. lati,long,area,rad,g,r,cpp)
304	!!! We check everything is OK.
305	PRINT *,"CHECK"
306	PRINT *,"--> mugaz = ",mugaz
307	PRINT *,"--> cpp = ",cpp
308	r = 8.314511E+0 * 1000.E+0 / mugaz
309	rcp = r / cpp
310
311
312	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
313	!!!! PLANETARY CONSTANTS !!!!
314	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
315
316	g = -99999.
317	PRINT *,'GRAVITY in m s-2 ?'
318	call getin("g",g)
319	IF (g.eq.-99999.) THEN
320	PRINT *,"STOP. I NEED g IN RUN.DEF."
321	STOP
322	ELSE
323	PRINT *,"--> g = ",g
324	ENDIF
325
326	!rad = -99999.
327	!PRINT *,'PLANETARY RADIUS in m ?'
328	!call getin("rad",rad)
329	!IF (rad.eq.-99999.) THEN
330	! PRINT *,"STOP. I NEED rad IN RUN.DEF."
331	! STOP
332	!ELSE
333	! PRINT *,"--> rad = ",rad
334	!ENDIF
335
336	daysec = -99999.
337	PRINT *,'LENGTH OF A DAY in s ?'
338	call getin("daysec",daysec)
339	IF (daysec.eq.-99999.) THEN
340	PRINT *,"STOP. I NEED daysec IN RUN.DEF."
341	STOP
342	ELSE
343	PRINT *,"--> daysec = ",daysec
344	ENDIF
345	omeg=4.*asin(1.)/(daysec)
346	PRINT *,"OK. FROM THIS I WORKED OUT:"
347	PRINT *,"--> omeg = ",omeg
348
349	year_day = -99999.
350	PRINT *,'LENGTH OF A YEAR in days ?'
351	call getin("year_day",year_day)
352	IF (year_day.eq.-99999.) THEN
353	PRINT *,"STOP. I NEED year_day IN RUN.DEF."
354	STOP
355	ELSE
356	PRINT *,"--> year_day = ",year_day
357	ENDIF
358
359	periastr = -99999.
360	PRINT *,'MIN DIST STAR-PLANET in AU ?'
361	call getin("periastr",periastr)
362	IF (periastr.eq.-99999.) THEN
363	PRINT *,"STOP. I NEED periastr IN RUN.DEF."
364	STOP
365	ELSE
366	PRINT *,"--> periastr = ",periastr
367	ENDIF
368
369	apoastr = -99999.
370	PRINT *,'MAX DIST STAR-PLANET in AU ?'
371	call getin("apoastr",apoastr)
372	IF (apoastr.eq.-99999.) THEN
373	PRINT *,"STOP. I NEED apoastr IN RUN.DEF."
374	STOP
375	ELSE
376	PRINT *,"--> apoastr = ",apoastr
377	ENDIF
378
379	peri_day = -99999.
380	PRINT *,'DATE OF PERIASTRON in days ?'
381	call getin("peri_day",peri_day)
382	IF (peri_day.eq.-99999.) THEN
383	PRINT *,"STOP. I NEED peri_day IN RUN.DEF."
384	STOP
385	ELSE IF (peri_day.gt.year_day) THEN
386	PRINT *,"STOP. peri_day.gt.year_day"
387	STOP
388	ELSE
389	PRINT *,"--> peri_day = ", peri_day
390	ENDIF
391
392	obliquit = -99999.
393	PRINT *,'OBLIQUITY in deg ?'
394	call getin("obliquit",obliquit)
395	IF (obliquit.eq.-99999.) THEN
396	PRINT *,"STOP. I NEED obliquit IN RUN.DEF."
397	STOP
398	ELSE
399	PRINT *,"--> obliquit = ",obliquit
400	ENDIF
401
402	psurf = -99999.
403	PRINT *,'SURFACE PRESSURE in Pa ?'
404	call getin("psurf",psurf)
405	IF (psurf.eq.-99999.) THEN
406	PRINT *,"STOP. I NEED psurf IN RUN.DEF."
407	STOP
408	ELSE
409	PRINT *,"--> psurf = ",psurf
410	ENDIF
411	!! we need reference pressures
412	pa=psurf/30.
413	preff=psurf ! these values are not needed in 1D (are you sure JL12)
414
415	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
416	!!!! END PLANETARY CONSTANTS !!!!
417	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
418
419	c Date et heure locale du debut du run
420	c ------------------------------------
421	c Date (en sols depuis le solstice de printemps) du debut du run
422	day0 = 0 ! default value for day0
423	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
424	call getin("day0",day0)
425	day=float(day0)
426	write(,) " day0 = ",day0
427	c Heure de demarrage
428	time=0 ! default value for time
429	write(,)'Initial local time (in hours, between 0 and 24)?'
430	call getin("time",time)
431	write(,)" time = ",time
432	time=time/24.E+0 ! convert time (hours) to fraction of sol
433
434
435	c Discretisation (Definition de la grille et des pas de temps)
436	c --------------
437	c
438	nlayer=nlayermx
439	nlevel=nlayer+1
440	nsoil=nsoilmx
441
442	day_step=48 ! default value for day_step
443	PRINT *,'Number of time steps per sol ?'
444	call getin("day_step",day_step)
445	write(,) " day_step = ",day_step
446
447
448	ecritphy=day_step ! default value for ecritphy
449	PRINT *,'Nunber of steps between writediagfi ?'
450	call getin("ecritphy",ecritphy)
451	write(,) " ecritphy = ",ecritphy
452
453	ndt=10 ! default value for ndt
454	PRINT *,'Number of sols to run ?'
455	call getin("ndt",ndt)
456	write(,) " ndt = ",ndt
457
458	ndt=ndt*day_step
459	dtphys=daysec/day_step
460	write(,)"-------------------------------------"
461	write(,)"-------------------------------------"
462	write(,)"--> Physical timestep is ",dtphys
463	write(,)"-------------------------------------"
464	write(,)"-------------------------------------"
465
466	c output spectrum?
467	write(,)"Output spectral OLR?"
468	specOLR=.false.
469	call getin("specOLR",specOLR)
470	write(,)" specOLR = ",specOLR
471
472	c
473	c pour le schema d'ondes de gravite
474	c ---------------------------------
475	zmea(1)=0.E+0
476	zstd(1)=0.E+0
477	zsig(1)=0.E+0
478	zgam(1)=0.E+0
479	zthe(1)=0.E+0
480
481	c Initialisation des traceurs
482	c ---------------------------
483
484	DO iq=1,nqmx
485	DO ilayer=1,nlayer
486	q(ilayer,iq) = 0.
487	ENDDO
488	ENDDO
489
490	DO iq=1,nqmx
491	qsurf(iq) = 0.
492	ENDDO
493
494	if (water) then
495	qzero1D = 0.0
496	qsurf(i_h2o_vap) = qzero1D
497	endif
498
499	c Initialisation pour prendre en compte les vents en 1-D
500	c ------------------------------------------------------
501	ptif=2.E+0omegsinlat(1)
502
503
504	c vent geostrophique
505	gru=10. ! default value for gru
506	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
507	call getin("u",gru)
508	write(,) " u = ",gru
509	grv=0. !default value for grv
510	PRINT *,'meridional northward component of the geostrophic',
511	&' wind (m/s) ?'
512	call getin("v",grv)
513	write(,) " v = ",grv
514
515	c Initialisation des vents au premier pas de temps
516	DO ilayer=1,nlayer
517	u(ilayer)=gru
518	v(ilayer)=grv
519	ENDDO
520
521	c energie cinetique turbulente
522	DO ilevel=1,nlevel
523	q2(ilevel)=0.E+0
524	ENDDO
525
526	c emissivity / surface co2 ice ( + h2o ice??)
527	c -------------------------------------------
528	emis=emissiv ! default value for emissivity
529	PRINT *,'Emissivity of bare ground ?'
530	call getin("emis",emis)
531	write(,) " emis = ",emis
532	emissiv=emis ! we do this so that condense_co2 sets things to the right
533	! value if there is no snow
534
535	if(i_co2_ice.gt.0)then
536	qsurf(i_co2_ice)=0 ! default value for co2ice
537	print*,'Initial CO2 ice on the surface (kg.m-2)'
538	call getin("co2ice",qsurf(i_co2_ice))
539	write(,) " co2ice = ",qsurf(i_co2_ice)
540	IF (qsurf(i_co2_ice).ge.1.E+0) THEN
541	! if we have some CO2 ice on the surface, change emissivity
542	if (lati(1).ge.0) then ! northern hemisphere
543	emis=emisice(1)
544	else ! southern hemisphere
545	emis=emisice(2)
546	endif
547	ENDIF
548	endif
549
550	c calcul des pressions et altitudes en utilisant les niveaux sigma
551	c ----------------------------------------------------------------
552
553	c Vertical Coordinates
554	c """"""""""""""""""""
555	hybrid=.true.
556	PRINT *,'Hybrid coordinates ?'
557	call getin("hybrid",hybrid)
558	write(,) " hybrid = ", hybrid
559
560
561	autozlevs=.false.
562	PRINT *,'Auto-discretise vertical levels ?'
563	call getin("autozlevs",autozlevs)
564	write(,) " autozlevs = ", autozlevs
565
566	pceil = psurf / 1000.0 ! Pascals
567	PRINT *,'Ceiling pressure (Pa) ?'
568	call getin("pceil",pceil)
569	write(,) " pceil = ", pceil
570
571	! Test of incompatibility:
572	! if autozlevs used, cannot have hybrid too
573
574	if (autozlevs.and.hybrid) then
575	print*,'Cannot use autozlevs and hybrid together!'
576	call abort
577	endif
578
579	if(autozlevs)then
580
581	open(91,file="z2sig.def",form='formatted')
582	read(91,*) Hscale
583	DO ilayer=1,nlayer-2
584	read(91,*) Hmax
585	enddo
586	close(91)
587
588
589	print*,'Hmax = ',Hmax,' km'
590	print*,'Auto-shifting Hscale to:'
591	! Hscale = Hmax / log(psurf/100.0)
592	Hscale = Hmax / log(psurf/pceil)
593	print*,'Hscale = ',Hscale,' km'
594
595	! none of this matters if we dont care about zlay
596
597	endif
598
599	call disvert
600
601	if(.not.autozlevs)then
602	! we want only the scale height from z2sig, in order to compute zlay
603	open(91,file="z2sig.def",form='formatted')
604	read(91,*) Hscale
605	close(91)
606	endif
607
608	! if(autozlevs)then
609	! open(94,file="Hscale.temp",form='formatted')
610	! read(94,*) Hscale
611	! close(94)
612	! endif
613
614	DO ilevel=1,nlevel
615	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
616	ENDDO
617
618	DO ilayer=1,nlayer
619	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
620	ENDDO
621
622
623
624	DO ilayer=1,nlayer
625	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
626	! & /g
627	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
628	ENDDO
629
630	! endif
631
632	c profil de temperature au premier appel
633	c --------------------------------------
634	pks=psurf**rcp
635
636	c altitude en km dans profile: on divise zlay par 1000
637	tmp1(0)=0.E+0
638	DO ilayer=1,nlayer
639	tmp1(ilayer)=zlay(ilayer)/1000.E+0
640	ENDDO
641	call profile(nlayer+1,tmp1,tmp2)
642
643	tsurf=tmp2(0)
644	DO ilayer=1,nlayer
645	temp(ilayer)=tmp2(ilayer)
646	ENDDO
647
648
649	c Initialisation albedo / inertie du sol
650	c --------------------------------------
651	albedodat(1)=0.2 ! default value for albedodat
652	PRINT *,'Albedo of bare ground ?'
653	call getin("albedo",albedodat(1))
654	write(,) " albedo = ",albedodat(1)
655
656	inertiedat(1,1)=400 ! default value for inertiedat
657	PRINT *,'Soil thermal inertia (SI) ?'
658	call getin("inertia",inertiedat(1,1))
659	write(,) " inertia = ",inertiedat(1,1)
660
661	! Initialize soil properties and temperature
662	! ------------------------------------------
663	volcapa=1.e6 ! volumetric heat capacity
664	DO isoil=1,nsoil
665	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
666	tsoil(isoil)=tsurf ! soil temperature
667	ENDDO
668
669	! Initialize depths
670	! -----------------
671	do isoil=0,nsoil-1
672	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
673	enddo
674	do isoil=1,nsoil
675	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
676	enddo
677
678
679	c Ecriture de "startfi"
680	c --------------------
681	c (Ce fichier sera aussitot relu au premier
682	c appel de "physiq", mais il est necessaire pour passer
683	c les variables purement physiques a "physiq"...
684
685	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
686	& dtphys,float(day0),
687	& time,tsurf,tsoil,emis,q2,qsurf,
688	& area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
689	& cloudfrac,totcloudfrac,hice)
690
691	c=======================================================================
692	c BOUCLE TEMPORELLE DU MODELE 1D
693	c=======================================================================
694
695	firstcall=.true.
696	lastcall=.false.
697
698	DO idt=1,ndt
699	IF (idt.eq.ndt) then !test
700	lastcall=.true.
701	call stellarlong(day*1.0,zls)
702	! write(103,) 'Ls=',zls180./pi
703	! write(103,) 'Lat=', lati(1)180./pi
704	! write(103,*) 'RunEnd - Atmos. Temp. File'
705	! write(103,*) 'RunEnd - Atmos. Temp. File'
706	! write(104,) 'Ls=',zls180./pi
707	! write(104,*) 'Lat=', lati(1)
708	! write(104,*) 'RunEnd - Atmos. Temp. File'
709	ENDIF
710
711	c calcul du geopotentiel
712	c ~~~~~~~~~~~~~~~~~~~~~
713
714
715	DO ilayer=1,nlayer
716
717	! if(autozlevs)then
718	! s(ilayer)=(play(ilayer)/psurf)**rcp
719	! else
720	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
721	! endif
722	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
723	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
724	ENDDO
725
726	! DO ilayer=1,nlayer
727	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
728	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
729	! ENDDO
730	phi(1)=pksh(1)(1.E+0-s(1))
731	DO ilayer=2,nlayer
732	phi(ilayer)=phi(ilayer-1)+
733	& pks(h(ilayer-1)+h(ilayer)).5E+0
734	& *(s(ilayer-1)-s(ilayer))
735
736	ENDDO
737
738	c appel de la physique
739	c --------------------
740
741
742	CALL physiq (1,llm,nqmx,
743	, firstcall,lastcall,
744	, day,time,dtphys,
745	, plev,play,phi,
746	, u, v,temp, q,
747	, w,
748	C - sorties
749	s du, dv, dtemp, dq,dpsurf,tracerdyn)
750
751
752	c evolution du vent : modele 1D
753	c -----------------------------
754
755	c la physique calcule les derivees temporelles de u et v.
756	c on y rajoute betement un effet Coriolis.
757	c
758	c DO ilayer=1,nlayer
759	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
760	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
761	c ENDDO
762
763	c Pour certain test : pas de coriolis a l'equateur
764	c if(lati(1).eq.0.) then
765	DO ilayer=1,nlayer
766	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
767	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
768	ENDDO
769	c end if
770	c
771	c
772	c Calcul du temps au pas de temps suivant
773	c ---------------------------------------
774	firstcall=.false.
775	time=time+dtphys/daysec
776	IF (time.gt.1.E+0) then
777	time=time-1.E+0
778	day=day+1
779	ENDIF
780
781	c calcul des vitesses et temperature au pas de temps suivant
782	c ----------------------------------------------------------
783
784	DO ilayer=1,nlayer
785	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
786	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
787	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
788	ENDDO
789
790	c calcul des pressions au pas de temps suivant
791	c ----------------------------------------------------------
792
793	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
794	DO ilevel=1,nlevel
795	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
796	ENDDO
797	DO ilayer=1,nlayer
798	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
799	ENDDO
800
801	c calcul traceur au pas de temps suivant
802	c --------------------------------------
803
804	DO iq = 1, nqmx
805	DO ilayer=1,nlayer
806	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
807	ENDDO
808	END DO
809
810	c ========================================================
811	c GESTION DES SORTIE
812	c ========================================================
813	if(saveprofile)then
814	OPEN(12,file='profile.out',form='formatted')
815	write(12,*) tsurf
816	DO ilayer=1,nlayermx
817	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
818	ENDDO
819	CLOSE(12)
820	endif
821
822
823	ENDDO ! fin de la boucle temporelle
824
825	c ========================================================
826	end !rcm1d
827
828	c***********************************************************************
829	c***********************************************************************
830	c Subroutines Bidons utilise seulement en 3D, mais
831	c necessaire a la compilation de rcm1d en 1D
832
833	subroutine gr_fi_dyn
834	RETURN
835	END
836
837	c***********************************************************************
838	c***********************************************************************
839
840	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format