Context Navigation

rcm1d.F @ 3567

Last change on this file since 3567 was 3562, checked in by mmaurice, 4 weeks ago

Generic PCM

1D restart operational: a restart1D.nc file is created that contains
psurf, tracers, winds and temperature profiles. A retartfi.nc file is
also created. Move those to and start1D.nc and startfi.nc and set
"restart" flag to .true. in rcm1d.def to restart from the files (also
make sure that day0 corresponds to the value in startfi.nc).

Property svn:executable set to *

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

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

Download in other formats:

Original Format