Context Navigation

rcm1d.F @ 3335

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

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

Download in other formats:

Original Format