Context Navigation

rcm1d.F @ 3580

Last change on this file since 3580 was 3576, checked in by jbclement, 6 months ago

COMMON:
Follow-up of r3574:

Small corrections to make it work with Git;
Addition of the functionality with the programs 'Generic newstart' and 'Generic start2archive';
Improvements of the visualization format;
Diplaying version control information of every sub-folder in the "trunk" instead of only the "trunk" (who can do more can do less).

JBC

Property svn:executable set to *

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

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

Download in other formats:

Original Format