Context Navigation

rcm1d.F @ 601

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