Context Navigation

rcm1d.F @ 1351

Last change on this file since 1351 was 1318, checked in by aslmd, 10 years ago
LMDZ.GENERIC. Two fixes for 1D version. First a conflict between comvert.h and planete_mod.F for ap bp preff (this one should be given a closer look, this redundancy is weird). Secondly output frequency was hardcoded to 1 (instead of ecritphy) in r1216.
Property svn:executable set to ``*
File size: 29.5 KB

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

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

Download in other formats:

Original Format