Context Navigation

rcm1d.F @ 3235

Last change on this file since 3235 was 3105, checked in by emillour, 2 years ago

Generic PCM:
Some corrections on the 1D models to keep up with recent code changes:

rcm1d: "noceanmx" no longer exists; slab size is always 2
kcm1d: naerkind is now a variable, so aerosol(:,:) needs be allocated.

Property svn:executable set to *

File size: 36.5 KB

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

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

Download in other formats:

Original Format