Context Navigation

rcm1d.F @ 1802

Last change on this file since 1802 was 1802, checked in by bclmd, 7 years ago
Adding photochemistry to LMDZ Generic
Property svn:executable set to ``*
File size: 33.1 KB

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

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

Download in other formats:

Original Format