Context Navigation

rcm1d.F @ 3523

Last change on this file since 3523 was 3515, checked in by emillour, 11 days ago
Generic PCM: Add a controle_descriptor to restartfi.nc file, describing contents of the controle array stored therein. Remove obsolete unused variables "lmixmin" and "emin_turb" in the process. EM
Property svn:executable set to ``*
File size: 36.6 KB

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

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

Download in other formats:

Original Format