Context Navigation

rcm1d.F @ 1243

Last change on this file since 1243 was 1216, checked in by emillour, 11 years ago

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

Property svn:executable set to *

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

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

Download in other formats:

Original Format