Context Navigation

rcm1d.F @ 724

Last change on this file since 724 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

Property svn:executable set to *

File size: 26.2 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 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	CALL inifis(1,llm,day0,daysec,dtphys,
302	. lati,long,area,rad,g,r,cpp)
303	!!! We check everything is OK.
304	PRINT *,"CHECK"
305	PRINT *,"--> mugaz = ",mugaz
306	PRINT *,"--> cpp = ",cpp
307	r = 8.314511E+0 * 1000.E+0 / mugaz
308	rcp = r / cpp
309
310
311	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312	!!!! PLANETARY CONSTANTS !!!!
313	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
314
315	g = -99999.
316	PRINT *,'GRAVITY in m s-2 ?'
317	call getin("g",g)
318	IF (g.eq.-99999.) THEN
319	PRINT *,"STOP. I NEED g IN RUN.DEF."
320	STOP
321	ELSE
322	PRINT *,"--> g = ",g
323	ENDIF
324
325	!rad = -99999.
326	!PRINT *,'PLANETARY RADIUS in m ?'
327	!call getin("rad",rad)
328	!IF (rad.eq.-99999.) THEN
329	! PRINT *,"STOP. I NEED rad IN RUN.DEF."
330	! STOP
331	!ELSE
332	! PRINT *,"--> rad = ",rad
333	!ENDIF
334
335	daysec = -99999.
336	PRINT *,'LENGTH OF A DAY in s ?'
337	call getin("daysec",daysec)
338	IF (daysec.eq.-99999.) THEN
339	PRINT *,"STOP. I NEED daysec IN RUN.DEF."
340	STOP
341	ELSE
342	PRINT *,"--> daysec = ",daysec
343	ENDIF
344	omeg=4.*asin(1.)/(daysec)
345	PRINT *,"OK. FROM THIS I WORKED OUT:"
346	PRINT *,"--> omeg = ",omeg
347
348	year_day = -99999.
349	PRINT *,'LENGTH OF A YEAR in days ?'
350	call getin("year_day",year_day)
351	IF (year_day.eq.-99999.) THEN
352	PRINT *,"STOP. I NEED year_day IN RUN.DEF."
353	STOP
354	ELSE
355	PRINT *,"--> year_day = ",year_day
356	ENDIF
357
358	periastr = -99999.
359	PRINT *,'MIN DIST STAR-PLANET in AU ?'
360	call getin("periastr",periastr)
361	IF (periastr.eq.-99999.) THEN
362	PRINT *,"STOP. I NEED periastr IN RUN.DEF."
363	STOP
364	ELSE
365	PRINT *,"--> periastr = ",periastr
366	ENDIF
367
368	apoastr = -99999.
369	PRINT *,'MAX DIST STAR-PLANET in AU ?'
370	call getin("apoastr",apoastr)
371	IF (apoastr.eq.-99999.) THEN
372	PRINT *,"STOP. I NEED apoastr IN RUN.DEF."
373	STOP
374	ELSE
375	PRINT *,"--> apoastr = ",apoastr
376	ENDIF
377
378	peri_day = -99999.
379	PRINT *,'DATE OF PERIASTRON in days ?'
380	call getin("peri_day",peri_day)
381	IF (peri_day.eq.-99999.) THEN
382	PRINT *,"STOP. I NEED peri_day IN RUN.DEF."
383	STOP
384	ELSE IF (peri_day.gt.year_day) THEN
385	PRINT *,"STOP. peri_day.gt.year_day"
386	STOP
387	ELSE
388	PRINT *,"--> peri_day = ", peri_day
389	ENDIF
390
391	obliquit = -99999.
392	PRINT *,'OBLIQUITY in deg ?'
393	call getin("obliquit",obliquit)
394	IF (obliquit.eq.-99999.) THEN
395	PRINT *,"STOP. I NEED obliquit IN RUN.DEF."
396	STOP
397	ELSE
398	PRINT *,"--> obliquit = ",obliquit
399	ENDIF
400
401	psurf = -99999.
402	PRINT *,'SURFACE PRESSURE in Pa ?'
403	call getin("psurf",psurf)
404	IF (psurf.eq.-99999.) THEN
405	PRINT *,"STOP. I NEED psurf IN RUN.DEF."
406	STOP
407	ELSE
408	PRINT *,"--> psurf = ",psurf
409	ENDIF
410	!! we need reference pressures
411	pa=psurf/30.
412	preff=psurf ! these values are not needed in 1D (are you sure JL12)
413
414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415	!!!! END PLANETARY CONSTANTS !!!!
416	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417
418	c Date et heure locale du debut du run
419	c ------------------------------------
420	c Date (en sols depuis le solstice de printemps) du debut du run
421	day0 = 0 ! default value for day0
422	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
423	call getin("day0",day0)
424	day=float(day0)
425	write(,) " day0 = ",day0
426	c Heure de demarrage
427	time=0 ! default value for time
428	write(,)'Initial local time (in hours, between 0 and 24)?'
429	call getin("time",time)
430	write(,)" time = ",time
431	time=time/24.E+0 ! convert time (hours) to fraction of sol
432
433
434	c Discretisation (Definition de la grille et des pas de temps)
435	c --------------
436	c
437	nlayer=nlayermx
438	nlevel=nlayer+1
439	nsoil=nsoilmx
440
441	day_step=48 ! default value for day_step
442	PRINT *,'Number of time steps per sol ?'
443	call getin("day_step",day_step)
444	write(,) " day_step = ",day_step
445
446
447	ecritphy=day_step ! default value for ecritphy
448	PRINT *,'Nunber of steps between writediagfi ?'
449	call getin("ecritphy",ecritphy)
450	write(,) " ecritphy = ",ecritphy
451
452	ndt=10 ! default value for ndt
453	PRINT *,'Number of sols to run ?'
454	call getin("ndt",ndt)
455	write(,) " ndt = ",ndt
456
457	ndt=ndt*day_step
458	dtphys=daysec/day_step
459	write(,)"-------------------------------------"
460	write(,)"-------------------------------------"
461	write(,)"--> Physical timestep is ",dtphys
462	write(,)"-------------------------------------"
463	write(,)"-------------------------------------"
464
465	c output spectrum?
466	write(,)"Output spectral OLR?"
467	specOLR=.false.
468	call getin("specOLR",specOLR)
469	write(,)" specOLR = ",specOLR
470
471	!!! AS12: what shall we do with this?
472	c Proprietes des poussiere aerosol
473	c --------------------------------
474	tauvis=0.2 ! default value for tauvis
475	print *,'Reference dust opacity at 700 Pa ?'
476	call getin("tauvis",tauvis)
477	write(,) " tauvis = ",tauvis
478
479	c
480	c pour le schema d'ondes de gravite
481	c ---------------------------------
482	zmea(1)=0.E+0
483	zstd(1)=0.E+0
484	zsig(1)=0.E+0
485	zgam(1)=0.E+0
486	zthe(1)=0.E+0
487
488	c Initialisation des traceurs
489	c ---------------------------
490
491	DO iq=1,nqmx
492	DO ilayer=1,nlayer
493	q(ilayer,iq) = 0.
494	ENDDO
495	ENDDO
496
497	DO iq=1,nqmx
498	qsurf(iq) = 0.
499	ENDDO
500
501	if (water) then
502	qzero1D = 0.0
503	qsurf(i_h2o_vap) = qzero1D
504	endif
505
506	c Initialisation pour prendre en compte les vents en 1-D
507	c ------------------------------------------------------
508	ptif=2.E+0omegsinlat(1)
509
510
511	c vent geostrophique
512	gru=10. ! default value for gru
513	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
514	call getin("u",gru)
515	write(,) " u = ",gru
516	grv=0. !default value for grv
517	PRINT *,'meridional northward component of the geostrophic',
518	&' wind (m/s) ?'
519	call getin("v",grv)
520	write(,) " v = ",grv
521
522	c Initialisation des vents au premier pas de temps
523	DO ilayer=1,nlayer
524	u(ilayer)=gru
525	v(ilayer)=grv
526	ENDDO
527
528	c energie cinetique turbulente
529	DO ilevel=1,nlevel
530	q2(ilevel)=0.E+0
531	ENDDO
532
533	c emissivity / surface co2 ice ( + h2o ice??)
534	c -------------------------------------------
535	emis=emissiv ! default value for emissivity
536	PRINT *,'Emissivity of bare ground ?'
537	call getin("emis",emis)
538	write(,) " emis = ",emis
539	emissiv=emis ! we do this so that condense_co2 sets things to the right
540	! value if there is no snow
541
542	if(i_co2_ice.gt.0)then
543	qsurf(i_co2_ice)=0 ! default value for co2ice
544	print*,'Initial CO2 ice on the surface (kg.m-2)'
545	call getin("co2ice",qsurf(i_co2_ice))
546	write(,) " co2ice = ",qsurf(i_co2_ice)
547	IF (qsurf(i_co2_ice).ge.1.E+0) THEN
548	! if we have some CO2 ice on the surface, change emissivity
549	if (lati(1).ge.0) then ! northern hemisphere
550	emis=emisice(1)
551	else ! southern hemisphere
552	emis=emisice(2)
553	endif
554	ENDIF
555	endif
556
557	c calcul des pressions et altitudes en utilisant les niveaux sigma
558	c ----------------------------------------------------------------
559
560	c Vertical Coordinates
561	c """"""""""""""""""""
562	hybrid=.true.
563	PRINT *,'Hybrid coordinates ?'
564	call getin("hybrid",hybrid)
565	write(,) " hybrid = ", hybrid
566
567
568	autozlevs=.false.
569	PRINT *,'Auto-discretise vertical levels ?'
570	call getin("autozlevs",autozlevs)
571	write(,) " autozlevs = ", autozlevs
572
573	pceil = psurf / 1000.0 ! Pascals
574	PRINT *,'Ceiling pressure (Pa) ?'
575	call getin("pceil",pceil)
576	write(,) " pceil = ", pceil
577
578	! Test of incompatibility:
579	! if autozlevs used, cannot have hybrid too
580
581	if (autozlevs.and.hybrid) then
582	print*,'Cannot use autozlevs and hybrid together!'
583	call abort
584	endif
585
586	if(autozlevs)then
587
588	open(91,file="z2sig.def",form='formatted')
589	read(91,*) Hscale
590	DO ilayer=1,nlayer-2
591	read(91,*) Hmax
592	enddo
593	close(91)
594
595
596	print*,'Hmax = ',Hmax,' km'
597	print*,'Auto-shifting Hscale to:'
598	! Hscale = Hmax / log(psurf/100.0)
599	Hscale = Hmax / log(psurf/pceil)
600	print*,'Hscale = ',Hscale,' km'
601
602	! none of this matters if we dont care about zlay
603
604	endif
605
606	call disvert
607
608	if(.not.autozlevs)then
609	! we want only the scale height from z2sig, in order to compute zlay
610	open(91,file="z2sig.def",form='formatted')
611	read(91,*) Hscale
612	close(91)
613	endif
614
615	! if(autozlevs)then
616	! open(94,file="Hscale.temp",form='formatted')
617	! read(94,*) Hscale
618	! close(94)
619	! endif
620
621	DO ilevel=1,nlevel
622	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
623	ENDDO
624
625	DO ilayer=1,nlayer
626	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
627	ENDDO
628
629
630
631	DO ilayer=1,nlayer
632	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
633	! & /g
634	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
635	ENDDO
636
637	! endif
638
639	c profil de temperature au premier appel
640	c --------------------------------------
641	pks=psurf**rcp
642
643	c altitude en km dans profile: on divise zlay par 1000
644	tmp1(0)=0.E+0
645	DO ilayer=1,nlayer
646	tmp1(ilayer)=zlay(ilayer)/1000.E+0
647	ENDDO
648	call profile(nlayer+1,tmp1,tmp2)
649
650	tsurf=tmp2(0)
651	DO ilayer=1,nlayer
652	temp(ilayer)=tmp2(ilayer)
653	ENDDO
654
655
656	c Initialisation albedo / inertie du sol
657	c --------------------------------------
658	albedodat(1)=0.2 ! default value for albedodat
659	PRINT *,'Albedo of bare ground ?'
660	call getin("albedo",albedodat(1))
661	write(,) " albedo = ",albedodat(1)
662
663	inertiedat(1,1)=400 ! default value for inertiedat
664	PRINT *,'Soil thermal inertia (SI) ?'
665	call getin("inertia",inertiedat(1,1))
666	write(,) " inertia = ",inertiedat(1,1)
667
668	! Initialize soil properties and temperature
669	! ------------------------------------------
670	volcapa=1.e6 ! volumetric heat capacity
671	DO isoil=1,nsoil
672	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
673	tsoil(isoil)=tsurf ! soil temperature
674	ENDDO
675
676	! Initialize depths
677	! -----------------
678	do isoil=0,nsoil-1
679	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
680	enddo
681	do isoil=1,nsoil
682	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
683	enddo
684
685
686	c Ecriture de "startfi"
687	c --------------------
688	c (Ce fichier sera aussitot relu au premier
689	c appel de "physiq", mais il est necessaire pour passer
690	c les variables purement physiques a "physiq"...
691
692	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
693	& dtphys,float(day0),
694	& time,tsurf,tsoil,emis,q2,qsurf,
695	& area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
696	& cloudfrac,totcloudfrac,hice)
697
698	c=======================================================================
699	c BOUCLE TEMPORELLE DU MODELE 1D
700	c=======================================================================
701
702	firstcall=.true.
703	lastcall=.false.
704
705	DO idt=1,ndt
706	IF (idt.eq.ndt) then !test
707	lastcall=.true.
708	call stellarlong(day*1.0,zls)
709	! write(103,) 'Ls=',zls180./pi
710	! write(103,) 'Lat=', lati(1)180./pi
711	! write(103,) 'Tau=', tauvis/700psurf
712	! write(103,*) 'RunEnd - Atmos. Temp. File'
713	! write(103,*) 'RunEnd - Atmos. Temp. File'
714	! write(104,) 'Ls=',zls180./pi
715	! write(104,*) 'Lat=', lati(1)
716	! write(104,) 'Tau=', tauvis/700psurf
717	! write(104,*) 'RunEnd - Atmos. Temp. File'
718	ENDIF
719
720	c calcul du geopotentiel
721	c ~~~~~~~~~~~~~~~~~~~~~
722
723
724	DO ilayer=1,nlayer
725
726	! if(autozlevs)then
727	! s(ilayer)=(play(ilayer)/psurf)**rcp
728	! else
729	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
730	! endif
731	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
732	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
733	ENDDO
734
735	! DO ilayer=1,nlayer
736	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
737	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
738	! ENDDO
739	phi(1)=pksh(1)(1.E+0-s(1))
740	DO ilayer=2,nlayer
741	phi(ilayer)=phi(ilayer-1)+
742	& pks(h(ilayer-1)+h(ilayer)).5E+0
743	& *(s(ilayer-1)-s(ilayer))
744
745	ENDDO
746
747	c appel de la physique
748	c --------------------
749
750
751	CALL physiq (1,llm,nqmx,
752	, firstcall,lastcall,
753	, day,time,dtphys,
754	, plev,play,phi,
755	, u, v,temp, q,
756	, w,
757	C - sorties
758	s du, dv, dtemp, dq,dpsurf,tracerdyn)
759
760
761	c evolution du vent : modele 1D
762	c -----------------------------
763
764	c la physique calcule les derivees temporelles de u et v.
765	c on y rajoute betement un effet Coriolis.
766	c
767	c DO ilayer=1,nlayer
768	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
769	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
770	c ENDDO
771
772	c Pour certain test : pas de coriolis a l'equateur
773	c if(lati(1).eq.0.) then
774	DO ilayer=1,nlayer
775	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
776	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
777	ENDDO
778	c end if
779	c
780	c
781	c Calcul du temps au pas de temps suivant
782	c ---------------------------------------
783	firstcall=.false.
784	time=time+dtphys/daysec
785	IF (time.gt.1.E+0) then
786	time=time-1.E+0
787	day=day+1
788	ENDIF
789
790	c calcul des vitesses et temperature au pas de temps suivant
791	c ----------------------------------------------------------
792
793	DO ilayer=1,nlayer
794	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
795	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
796	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
797	ENDDO
798
799	c calcul des pressions au pas de temps suivant
800	c ----------------------------------------------------------
801
802	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
803	DO ilevel=1,nlevel
804	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
805	ENDDO
806	DO ilayer=1,nlayer
807	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
808	ENDDO
809
810	c calcul traceur au pas de temps suivant
811	c --------------------------------------
812
813	DO iq = 1, nqmx
814	DO ilayer=1,nlayer
815	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
816	ENDDO
817	END DO
818
819	c ========================================================
820	c GESTION DES SORTIE
821	c ========================================================
822	if(saveprofile)then
823	OPEN(12,file='profile.out',form='formatted')
824	write(12,*) tsurf
825	DO ilayer=1,nlayermx
826	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
827	ENDDO
828	CLOSE(12)
829	endif
830
831
832	ENDDO ! fin de la boucle temporelle
833
834	c ========================================================
835	end !rcm1d
836
837	c***********************************************************************
838	c***********************************************************************
839	c Subroutines Bidons utilise seulement en 3D, mais
840	c necessaire a la compilation de rcm1d en 1D
841
842	subroutine gr_fi_dyn
843	RETURN
844	END
845
846	c***********************************************************************
847	c***********************************************************************
848
849	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format