Context Navigation

rcm1d.F @ 4027

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

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

Download in other formats:

Original Format