Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rcm1d.F.delete @ 264

Last change on this file since 264 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

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

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

Download in other formats:

Original Format