Context Navigation

rcm1d.F @ 3197

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

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

Download in other formats:

Original Format