Context Navigation

rcm1d.F @ 1145

Last change on this file since 1145 was 970, checked in by aslmd, 12 years ago
LMZ.GENERIC rcm1d added a warning about tracer number inconsistency to avoid problems.
Property svn:executable set to ``*
File size: 27.6 KB

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

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

Download in other formats:

Original Format