Context Navigation

rcm1d.F @ 1990

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

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

Download in other formats:

Original Format