Context Navigation

rcm1d.F @ 2908

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

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

Download in other formats:

Original Format