Context Navigation

rcm1d.F @ 3000

Last change on this file since 3000 was 2972, checked in by emillour, 3 years ago
Generic PCM: Make number of scatterers fully dynamic (i.e. set in callphys.def and no longer by compilation option "-s #"). One should now specify naerkind = # in callphys.def (default is 0). EM
Property svn:executable set to ``*
File size: 36.5 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 surfdat_h, only: albedodat, phisfi, dryness,
10	& zmea, zstd, zsig, zgam, zthe,
11	& emissiv, emisice, iceradius,
12	& dtemisice
13	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
14	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
15	use phyredem, only: physdem0,physdem1
16	use geometry_mod, only: init_geometry
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	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) ! surface layer
98	REAL q2(llm+1) ! Turbulent Kinetic Energy
99	REAL zlay(llm) ! altitude estimee dans les couches (km)
100
101	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
102	REAL du(llm),dv(llm),dtemp(llm)
103	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
104	REAL dpsurf(1)
105	REAL,ALLOCATABLE :: dq(:,:)
106	REAL,ALLOCATABLE :: dqdyn(:,:)
107
108	c Various intermediate variables
109	! INTEGER thermo
110	REAL zls
111	REAL phi(llm),h(llm),s(llm)
112	REAL pks, ptif, w(llm)
113	INTEGER ierr, aslun
114	REAL tmp1(0:llm),tmp2(0:llm)
115	integer :: nq !=1 ! number of tracers
116
117	character*2 str2
118	character (len=7) :: str7
119	character(len=44) :: txt
120	character(len=44) :: tracer_profile_file_name
121
122	logical oldcompare, earthhack,saveprofile
123
124	! added by RW for zlay computation
125	real Hscale, Hmax, rho, dz
126
127	! added by RW for autozlevs computation
128	logical autozlevs
129	real nu, xx, pMIN, zlev, Htop
130	real logplevs(llm)
131
132	! added by BC
133	REAL cloudfrac(1,llm)
134	REAL hice(1),totcloudfrac(1)
135
136	! added by BC for ocean
137	real rnat(1)
138	REAL tslab(1,noceanmx),tsea_ice(1),sea_ice(1)
139	real pctsrf_sic(1)
140
141	! added by BC to accelerate convergence
142	logical accelerate_temperature
143	real factor_acceleration
144
145	! added by AS to avoid the use of adv trac common
146	character*30,allocatable :: nametmp(:) !
147
148	real :: latitude(1), longitude(1), cell_area(1)
149
150	! added by JVO and YJ to read modern traceur.def
151	character(len=500) :: line ! to store a line of text
152	LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def
153
154	c=======================================================================
155	c INITIALISATION
156	c=======================================================================
157	! check if 'rcm1d.def' file is around
158	open(90,file='rcm1d.def',status='old',form='formatted',
159	& iostat=ierr)
160	if (ierr.ne.0) then
161	write(,) 'Cannot find required file "rcm1d.def"'
162	write(,) 'which should contain some input parameters'
163	write(,) ' ... might as well stop here ...'
164	stop
165	else
166	close(90)
167	endif
168
169	! now, run.def is needed anyway. so we create a dummy temporary one
170	! for ioipsl to work. if a run.def is already here, stop the
171	! program and ask the user to do a bit of cleaning
172	open(90,file='run.def',status='old',form='formatted',
173	& iostat=ierr)
174	if (ierr.eq.0) then
175	close(90)
176	write(,) 'There is already a run.def file.'
177	write(,) 'This is not compatible with 1D runs.'
178	write(,) 'Please remove the file and restart the run.'
179	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
180	stop
181	else
182	call system('touch run.def')
183	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
184	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
185	endif
186
187	! read nq from traceur.def
188	open(90,file='traceur.def',status='old',form='formatted',
189	& iostat=ierr)
190	if (ierr.eq.0) then
191	! read number of tracers:
192	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
193	READ(90,'(A)') line
194	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
195	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
196	ELSE
197	moderntracdef = .true.
198	DO
199	READ(90,'(A)',iostat=ierr) line
200	IF (ierr.eq.0) THEN
201	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
202	READ(line,*,iostat=ierr) nq
203	EXIT
204	ENDIF
205	ELSE ! If pb, or if reached EOF without having found nbtr
206	write(,) "rcm1d: error reading number of tracers"
207	write(,) " (first line of traceur.def) "
208	stop
209	ENDIF
210	ENDDO
211	ENDIF ! if modern or standard traceur.def
212	else
213	nq=0
214	endif
215	close(90)
216
217	! Initialize dimphy module
218	call init_dimphy(1,llm)
219
220	! now initialize arrays using phys_state_var_init
221	! but first initialise naerkind (from callphys.def)
222	naerkind=0 !default
223	call getin("naerkind",naerkind)
224
225	call phys_state_var_init(nq)
226
227	saveprofile=.false.
228	saveprofile=.true.
229
230	c ----------------------------------------
231	c Default values (corresponding to Mars)
232	c ----------------------------------------
233
234	pi=2.E+0*asin(1.E+0)
235
236	c Parametres Couche limite et Turbulence
237	c --------------------------------------
238	z0 = 1.e-2 ! surface roughness (m) ~0.01
239	emin_turb = 1.e-6 ! energie minimale ~1.e-8
240	lmixmin = 30 ! longueur de melange ~100
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
887	inertiedat(1,1)=400 ! default value for inertiedat
888	PRINT *,'Soil thermal inertia (SI) ?'
889	call getin("inertia",inertiedat(1,1))
890	write(,) " inertia = ",inertiedat(1,1)
891
892	! Initialize soil properties and temperature
893	! ------------------------------------------
894	volcapa=1.e6 ! volumetric heat capacity
895	DO isoil=1,nsoil
896	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
897	tsoil(isoil)=tsurf(1) ! soil temperature
898	ENDDO
899
900	! Initialize depths
901	! -----------------
902	do isoil=0,nsoil-1
903	mlayer(isoil)=3.e-2(2.*(isoil-0.5)) ! mid-layer depth
904	enddo
905	do isoil=1,nsoil
906	layer(isoil)=3.e-2(2.*(isoil-1)) ! layer depth
907	enddo
908
909	! Initialize cloud fraction and oceanic ice
910	! -----------------------------------------
911	hice=0.
912	do ilayer=1,llm
913	cloudfrac(1,ilayer)=0.
914	enddo
915	totcloudfrac=0.
916
917	! Initialize slab ocean
918	! -----------------
919	rnat=1. ! default value for rnat
920	if(inertiedat(1,1).GE.10000.)then
921	rnat=0.
922	endif
923	if(ok_slab_ocean)then
924	rnat=0.
925	tslab(1,1)=tsurf(1)
926	tslab(1,2)=tsurf(1)
927	tsea_ice=tsurf
928	pctsrf_sic=0.
929	sea_ice=0.
930	endif
931
932
933	! Initialize chemical species
934	! -----------------
935	#ifndef MESOSCALE
936	if(tracer.and.photochem) then
937	call initracer(1,nq)
938	allocate(nametmp(nq))
939	nametmp(1:nq)=tname(1:nq)
940	call inichim_1D(nlayer, nq, q, qsurf, play, 0, 0)
941	tname(1:nq)=nametmp(1:nq)
942	noms(1:nq)=nametmp(1:nq)
943	endif ! tracer and photochem
944	#endif
945
946
947	c Write a "startfi" file
948	c --------------------
949	c This file will be read during the first call to "physiq".
950	c It is needed to transfert physics variables to "physiq"...
951
952	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
953	& dtphys,real(day0),time,cell_area,
954	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
955	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
956	& dtphys,time,
957	& tsurf,tsoil,emis,q2,qsurf,
958	& cloudfrac,totcloudfrac,hice,
959	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
960
961	c=======================================================================
962	c BOUCLE TEMPORELLE DU MODELE 1D
963	c=======================================================================
964
965	firstcall=.true.
966	lastcall=.false.
967
968	DO idt=1,ndt
969	IF (idt.eq.ndt) then !test
970	lastcall=.true.
971	call stellarlong(day*1.0,zls)
972	! write(103,) 'Ls=',zls180./pi
973	! write(103,) 'Lat=', latitude(1)180./pi
974	! write(103,*) 'RunEnd - Atmos. Temp. File'
975	! write(103,*) 'RunEnd - Atmos. Temp. File'
976	! write(104,) 'Ls=',zls180./pi
977	! write(104,*) 'Lat=', latitude(1)
978	! write(104,*) 'RunEnd - Atmos. Temp. File'
979	ENDIF
980
981	c calcul du geopotentiel
982	c ~~~~~~~~~~~~~~~~~~~~~
983
984
985	DO ilayer=1,nlayer
986
987	! if(autozlevs)then
988	! s(ilayer)=(play(ilayer)/psurf)**rcp
989	! else
990	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
991	! endif
992	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
993	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
994	ENDDO
995
996	! DO ilayer=1,nlayer
997	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
998	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
999	! ENDDO
1000	phi(1)=pksh(1)(1.E+0-s(1))
1001	DO ilayer=2,nlayer
1002	phi(ilayer)=phi(ilayer-1)+
1003	& pks(h(ilayer-1)+h(ilayer)).5E+0
1004	& *(s(ilayer-1)-s(ilayer))
1005
1006	ENDDO
1007
1008	c appel de la physique
1009	c --------------------
1010
1011
1012	CALL physiq (1,llm,nq,
1013	, firstcall,lastcall,
1014	, day,time,dtphys,
1015	, plev,play,phi,
1016	, u, v,temp, q,
1017	, w,
1018	C - sorties
1019	s du, dv, dtemp, dq,dpsurf)
1020
1021
1022	c evolution du vent : modele 1D
1023	c -----------------------------
1024
1025	c la physique calcule les derivees temporelles de u et v.
1026	c on y rajoute betement un effet Coriolis.
1027	c
1028	c DO ilayer=1,nlayer
1029	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
1030	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
1031	c ENDDO
1032
1033	c Pour certain test : pas de coriolis a l'equateur
1034	c if(latitude(1).eq.0.) then
1035	DO ilayer=1,nlayer
1036	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
1037	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
1038	ENDDO
1039	c end if
1040	c
1041	c
1042	c Calcul du temps au pas de temps suivant
1043	c ---------------------------------------
1044	firstcall=.false.
1045	time=time+dtphys/daysec
1046	IF (time.gt.1.E+0) then
1047	time=time-1.E+0
1048	day=day+1
1049	ENDIF
1050
1051	c calcul des vitesses et temperature au pas de temps suivant
1052	c ----------------------------------------------------------
1053
1054	DO ilayer=1,nlayer
1055	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
1056	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
1057	if(accelerate_temperature) then
1058	temp(ilayer)=temp(ilayer)+dtphysdtemp(ilayer)
1059	& max(1.0,play(ilayer)1e-5)*factor_acceleration
1060	else
1061	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
1062	endif
1063	ENDDO
1064
1065	c calcul des pressions au pas de temps suivant
1066	c ----------------------------------------------------------
1067
1068	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
1069	DO ilevel=1,nlevel
1070	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
1071	ENDDO
1072	DO ilayer=1,nlayer
1073	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
1074	ENDDO
1075
1076	c calcul traceur au pas de temps suivant
1077	c --------------------------------------
1078
1079	DO iq = 1, nq
1080	DO ilayer=1,nlayer
1081	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
1082	ENDDO
1083	END DO
1084
1085	c ========================================================
1086	c GESTION DES SORTIE
1087	c ========================================================
1088	if(saveprofile)then
1089	OPEN(12,file='profile.out',form='formatted')
1090	write(12,*) tsurf
1091	DO ilayer=1,llm
1092	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
1093	ENDDO
1094	CLOSE(12)
1095	endif
1096
1097
1098	ENDDO ! fin de la boucle temporelle
1099
1100	write(,) "rcm1d: Everything is cool."
1101
1102	c ========================================================
1103	end !rcm1d
1104
1105

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

Download in other formats:

Original Format