Context Navigation

rcm1d.F @ 2621

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

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

Download in other formats:

Original Format