Context Navigation

rcm1d.F @ 2276

Last change on this file since 2276 was 2076, checked in by jvatant, 6 years ago
Fix a pb of dummy init of r and cpp in rcm1d, leading to NaNs? with recent picky ifort compilers --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=-9999. ! dummy init for inifis, will be rewrite later on
520	r=-9999. ! dummy init for inifis, will be rewrite later on
521	CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
522	. latitude,longitude,cell_area,rad,g,r,cpp)
523
524	nsoil=nsoilmx
525	allocate(tsoil(nsoilmx))
526	!! those are defined in comsoil_h.F90
527	IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
528	IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
529	IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
530
531	! At this point, both getin() and getin_p() functions have been used,
532	! and the run.def file can be removed.
533	call system("rm -rf run.def")
534
535	!!! We check everything is OK.
536	PRINT *,"CHECK"
537	PRINT *,"--> mugaz = ",mugaz
538	PRINT *,"--> cpp = ",cpp
539	r = 8.314511E+0 * 1000.E+0 / mugaz
540	rcp = r / cpp
541
542	c output spectrum?
543	write(,)"Output spectral OLR?"
544	specOLR=.false.
545	call getin("specOLR",specOLR)
546	write(,)" specOLR = ",specOLR
547
548	c
549	c pour le schema d'ondes de gravite
550	c ---------------------------------
551	zmea(1)=0.E+0
552	zstd(1)=0.E+0
553	zsig(1)=0.E+0
554	zgam(1)=0.E+0
555	zthe(1)=0.E+0
556
557	c Initialisation des traceurs
558	c ---------------------------
559
560	DO iq=1,nq
561	DO ilayer=1,nlayer
562	q(ilayer,iq) = 0.
563	ENDDO
564	ENDDO
565
566	DO iq=1,nq
567	qsurf(iq) = 0.
568	ENDDO
569
570	if (tracer) then ! Initialize tracers here.
571
572	write(,) "rcm1d : initializing tracers profiles"
573
574	do iq=1,nq
575
576	txt=""
577	write(txt,"(a)") tname(iq)
578	write(,)" tracer:",trim(txt)
579
580	! CO2
581	if (txt.eq."co2_ice") then
582	q(:,iq)=0. ! kg/kg of atmosphere
583	qsurf(iq)=0. ! kg/m2 at the surface
584	! Look for a "profile_co2_ice" input file
585	open(91,file='profile_co2_ice',status='old',
586	& form='formatted',iostat=ierr)
587	if (ierr.eq.0) then
588	read(91,*) qsurf(iq)
589	do ilayer=1,nlayer
590	read(91,*) q(ilayer,iq)
591	enddo
592	else
593	write(,) "No profile_co2_ice file!"
594	endif
595	close(91)
596	endif ! of if (txt.eq."co2")
597
598	! WATER VAPOUR
599	if (txt.eq."h2o_vap") then
600	q(:,iq)=0. ! kg/kg of atmosphere
601	qsurf(iq)=0. ! kg/m2 at the surface
602	! Look for a "profile_h2o_vap" input file
603	open(91,file='profile_h2o_vap',status='old',
604	& form='formatted',iostat=ierr)
605	if (ierr.eq.0) then
606	read(91,*) qsurf(iq)
607	do ilayer=1,nlayer
608	read(91,*) q(ilayer,iq)
609	enddo
610	else
611	write(,) "No profile_h2o_vap file!"
612	endif
613	close(91)
614	endif ! of if (txt.eq."h2o_vap")
615
616	! WATER ICE
617	if (txt.eq."h2o_ice") then
618	q(:,iq)=0. ! kg/kg of atmosphere
619	qsurf(iq)=0. ! kg/m2 at the surface
620	! Look for a "profile_h2o_ice" input file
621	open(91,file='profile_h2o_ice',status='old',
622	& form='formatted',iostat=ierr)
623	if (ierr.eq.0) then
624	read(91,*) qsurf(iq)
625	do ilayer=1,nlayer
626	read(91,*) q(ilayer,iq)
627	enddo
628	else
629	write(,) "No profile_h2o_ice file!"
630	endif
631	close(91)
632	endif ! of if (txt.eq."h2o_ice")
633
634	enddo ! of do iq=1,nq
635
636	endif ! of tracer
637
638	c Initialisation pour prendre en compte les vents en 1-D
639	c ------------------------------------------------------
640	ptif=2.E+0omegsinlat(1)
641
642
643	c vent geostrophique
644	gru=10. ! default value for gru
645	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
646	call getin("u",gru)
647	write(,) " u = ",gru
648	grv=0. !default value for grv
649	PRINT *,'meridional northward component of the geostrophic',
650	&' wind (m/s) ?'
651	call getin("v",grv)
652	write(,) " v = ",grv
653
654	! To be clean, also set vertical winds to zero
655	w(1:nlayer)=0
656
657	c Initialisation des vents au premier pas de temps
658	DO ilayer=1,nlayer
659	u(ilayer)=gru
660	v(ilayer)=grv
661	ENDDO
662
663	c energie cinetique turbulente
664	DO ilevel=1,nlevel
665	q2(ilevel)=0.E+0
666	ENDDO
667
668	c emissivity / surface co2 ice ( + h2o ice??)
669	c -------------------------------------------
670	emis(1)=emissiv ! default value for emissivity
671	PRINT *,'Emissivity of bare ground ?'
672	call getin("emis",emis(1))
673	write(,) " emis = ",emis(1)
674	emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
675	! value if there is no snow
676
677	if(i_co2_ice.gt.0)then
678	qsurf(i_co2_ice)=0 ! default value for co2ice
679	print*,'Initial CO2 ice on the surface (kg.m-2)'
680	call getin("co2ice",qsurf(i_co2_ice))
681	write(,) " co2ice = ",qsurf(i_co2_ice)
682	IF (qsurf(i_co2_ice).ge.1.E+0) THEN
683	! if we have some CO2 ice on the surface, change emissivity
684	if (latitude(1).ge.0) then ! northern hemisphere
685	emis(1)=emisice(1)
686	else ! southern hemisphere
687	emis(1)=emisice(2)
688	endif
689	ENDIF
690	endif
691
692	c calcul des pressions et altitudes en utilisant les niveaux sigma
693	c ----------------------------------------------------------------
694
695	c Vertical Coordinates
696	c """"""""""""""""""""
697	hybrid=.true.
698	PRINT *,'Hybrid coordinates ?'
699	call getin("hybrid",hybrid)
700	write(,) " hybrid = ", hybrid
701
702
703	autozlevs=.false.
704	PRINT *,'Auto-discretise vertical levels ?'
705	call getin("autozlevs",autozlevs)
706	write(,) " autozlevs = ", autozlevs
707
708	pceil = psurf / 1000.0 ! Pascals
709	PRINT *,'Ceiling pressure (Pa) ?'
710	call getin("pceil",pceil)
711	write(,) " pceil = ", pceil
712
713	! Test of incompatibility:
714	! if autozlevs used, cannot have hybrid too
715
716	if (autozlevs.and.hybrid) then
717	print*,'Cannot use autozlevs and hybrid together!'
718	call abort
719	endif
720
721	if(autozlevs)then
722
723	open(91,file="z2sig.def",form='formatted')
724	read(91,*) Hscale
725	DO ilayer=1,nlayer-2
726	read(91,*) Hmax
727	enddo
728	close(91)
729
730
731	print*,'Hmax = ',Hmax,' km'
732	print*,'Auto-shifting Hscale to:'
733	! Hscale = Hmax / log(psurf/100.0)
734	Hscale = Hmax / log(psurf/pceil)
735	print*,'Hscale = ',Hscale,' km'
736
737	! none of this matters if we dont care about zlay
738
739	endif
740
741	call disvert
742	! now that disvert has been called, initialize module vertical_layers_mod
743	call init_vertical_layers(nlayer,preff,scaleheight,
744	& ap,bp,aps,bps,presnivs,pseudoalt)
745
746	if(.not.autozlevs)then
747	! we want only the scale height from z2sig, in order to compute zlay
748	open(91,file="z2sig.def",form='formatted')
749	read(91,*) Hscale
750	close(91)
751	endif
752
753	! if(autozlevs)then
754	! open(94,file="Hscale.temp",form='formatted')
755	! read(94,*) Hscale
756	! close(94)
757	! endif
758
759	DO ilevel=1,nlevel
760	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
761	ENDDO
762
763	DO ilayer=1,nlayer
764	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
765	ENDDO
766
767
768
769	DO ilayer=1,nlayer
770	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
771	! & /g
772	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
773	ENDDO
774
775	! endif
776
777	c profil de temperature au premier appel
778	c --------------------------------------
779	pks=psurf**rcp
780
781	c altitude en km dans profile: on divise zlay par 1000
782	tmp1(0)=0.E+0
783	DO ilayer=1,nlayer
784	tmp1(ilayer)=zlay(ilayer)/1000.E+0
785	ENDDO
786	call profile(nlayer+1,tmp1,tmp2)
787
788	tsurf(1)=tmp2(0)
789	DO ilayer=1,nlayer
790	temp(ilayer)=tmp2(ilayer)
791	ENDDO
792	print*,"check"
793	PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
794	PRINT*,"INPUT TEMPERATURE PROFILE",temp
795
796	c Initialisation albedo / inertie du sol
797	c --------------------------------------
798	albedodat(1)=0.2 ! default value for albedodat
799	PRINT *,'Albedo of bare ground ?'
800	call getin("albedo",albedodat(1))
801	write(,) " albedo = ",albedodat(1)
802
803	inertiedat(1,1)=400 ! default value for inertiedat
804	PRINT *,'Soil thermal inertia (SI) ?'
805	call getin("inertia",inertiedat(1,1))
806	write(,) " inertia = ",inertiedat(1,1)
807
808	! Initialize soil properties and temperature
809	! ------------------------------------------
810	volcapa=1.e6 ! volumetric heat capacity
811	DO isoil=1,nsoil
812	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
813	tsoil(isoil)=tsurf(1) ! soil temperature
814	ENDDO
815
816	! Initialize depths
817	! -----------------
818	do isoil=0,nsoil-1
819	mlayer(isoil)=3.e-2(2.*(isoil-0.5)) ! mid-layer depth
820	enddo
821	do isoil=1,nsoil
822	layer(isoil)=3.e-2(2.*(isoil-1)) ! layer depth
823	enddo
824
825	! Initialize cloud fraction and oceanic ice
826	! -----------------------------------------
827	hice=0.
828	do ilayer=1,llm
829	cloudfrac(1,ilayer)=0.
830	enddo
831	totcloudfrac=0.
832
833	! Initialize slab ocean
834	! -----------------
835	rnat=1. ! default value for rnat
836	if(inertiedat(1,1).GE.10000.)then
837	rnat=0.
838	endif
839	if(ok_slab_ocean)then
840	rnat=0.
841	tslab(1,1)=tsurf(1)
842	tslab(1,2)=tsurf(1)
843	tsea_ice=tsurf
844	pctsrf_sic=0.
845	sea_ice=0.
846	endif
847
848
849	! Initialize chemical species
850	! -----------------
851	#ifndef MESOSCALE
852	if(tracer.and.photochem) then
853	call initracer(1,nq,tname)
854	allocate(nametmp(nq))
855	nametmp(1:nq)=tname(1:nq)
856	call inichim_1D(nq, q, qsurf, psurf, 0, 0)
857	tname(1:nq)=nametmp(1:nq)
858	noms(1:nq)=nametmp(1:nq)
859	endif ! tracer and photochem
860	#endif
861
862
863	c Write a "startfi" file
864	c --------------------
865	c This file will be read during the first call to "physiq".
866	c It is needed to transfert physics variables to "physiq"...
867
868	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
869	& dtphys,real(day0),time,cell_area,
870	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
871	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
872	& dtphys,time,
873	& tsurf,tsoil,emis,q2,qsurf,
874	& cloudfrac,totcloudfrac,hice,
875	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
876
877	c=======================================================================
878	c BOUCLE TEMPORELLE DU MODELE 1D
879	c=======================================================================
880
881	firstcall=.true.
882	lastcall=.false.
883
884	DO idt=1,ndt
885	IF (idt.eq.ndt) then !test
886	lastcall=.true.
887	call stellarlong(day*1.0,zls)
888	! write(103,) 'Ls=',zls180./pi
889	! write(103,) 'Lat=', latitude(1)180./pi
890	! write(103,*) 'RunEnd - Atmos. Temp. File'
891	! write(103,*) 'RunEnd - Atmos. Temp. File'
892	! write(104,) 'Ls=',zls180./pi
893	! write(104,*) 'Lat=', latitude(1)
894	! write(104,*) 'RunEnd - Atmos. Temp. File'
895	ENDIF
896
897	c calcul du geopotentiel
898	c ~~~~~~~~~~~~~~~~~~~~~
899
900
901	DO ilayer=1,nlayer
902
903	! if(autozlevs)then
904	! s(ilayer)=(play(ilayer)/psurf)**rcp
905	! else
906	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
907	! endif
908	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
909	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
910	ENDDO
911
912	! DO ilayer=1,nlayer
913	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
914	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
915	! ENDDO
916	phi(1)=pksh(1)(1.E+0-s(1))
917	DO ilayer=2,nlayer
918	phi(ilayer)=phi(ilayer-1)+
919	& pks(h(ilayer-1)+h(ilayer)).5E+0
920	& *(s(ilayer-1)-s(ilayer))
921
922	ENDDO
923
924	c appel de la physique
925	c --------------------
926
927
928	CALL physiq (1,llm,nq,
929	. tname,
930	, firstcall,lastcall,
931	, day,time,dtphys,
932	, plev,play,phi,
933	, u, v,temp, q,
934	, w,
935	C - sorties
936	s du, dv, dtemp, dq,dpsurf)
937
938
939	c evolution du vent : modele 1D
940	c -----------------------------
941
942	c la physique calcule les derivees temporelles de u et v.
943	c on y rajoute betement un effet Coriolis.
944	c
945	c DO ilayer=1,nlayer
946	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
947	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
948	c ENDDO
949
950	c Pour certain test : pas de coriolis a l'equateur
951	c if(latitude(1).eq.0.) then
952	DO ilayer=1,nlayer
953	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
954	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
955	ENDDO
956	c end if
957	c
958	c
959	c Calcul du temps au pas de temps suivant
960	c ---------------------------------------
961	firstcall=.false.
962	time=time+dtphys/daysec
963	IF (time.gt.1.E+0) then
964	time=time-1.E+0
965	day=day+1
966	ENDIF
967
968	c calcul des vitesses et temperature au pas de temps suivant
969	c ----------------------------------------------------------
970
971	DO ilayer=1,nlayer
972	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
973	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
974	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
975	ENDDO
976
977	c calcul des pressions au pas de temps suivant
978	c ----------------------------------------------------------
979
980	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
981	DO ilevel=1,nlevel
982	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
983	ENDDO
984	DO ilayer=1,nlayer
985	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
986	ENDDO
987
988	c calcul traceur au pas de temps suivant
989	c --------------------------------------
990
991	DO iq = 1, nq
992	DO ilayer=1,nlayer
993	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
994	ENDDO
995	END DO
996
997	c ========================================================
998	c GESTION DES SORTIE
999	c ========================================================
1000	if(saveprofile)then
1001	OPEN(12,file='profile.out',form='formatted')
1002	write(12,*) tsurf
1003	DO ilayer=1,llm
1004	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
1005	ENDDO
1006	CLOSE(12)
1007	endif
1008
1009
1010	ENDDO ! fin de la boucle temporelle
1011
1012	write(,) "rcm1d: Everything is cool."
1013
1014	c ========================================================
1015	end !rcm1d
1016
1017

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

Download in other formats:

Original Format