Context Navigation

rcm1d.F @ 3574

Last change on this file since 3574 was 3574, checked in by jbclement, 2 weeks ago
COMMON: The compilation generates a Fortran subroutine to track compilation and version (SVN or Git) details through the executable file. Put the argument 'version' as an option when executing the code to display these information and create a file "version_diff.txt" containing the diff result. It can work with every program but it has been implemented only for: 'gcm', 'parallel gcm', 'pem', '1D Mars PCM', 'Mars newstart', 'Mars start2archive' and '1D Generic PCM'. JBC
Property svn:executable set to ``*
File size: 44.1 KB

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

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

Download in other formats:

Original Format