Context Navigation

rcm1d.F @ 600

Last change on this file since 600 was 591, checked in by aslmd, 13 years ago
LMDZ.GENERIC: sample files for 1D jupiter-like calculations in deftank. cleaned and compliant withy latest modifications.
Property svn:executable set to ``*
File size: 25.8 KB

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

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

Download in other formats:

Original Format