Context Navigation

scm.f90 @ 5623

Last change on this file since 5623 was 5609, checked in by aborella, 3 months ago

changed treatment of prognostic variables for prognostic clouds
adapted sedimentation and autoconversion for prognostic cirrus clouds
cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
simplified assumptions for ice crystals deposition / sublimation
first version of the coupling between prognostic cirrus clouds and deep convection
added persistent contrail cirrus clouds in radiative diagnostics

File size: 40.5 KB

Line
1	SUBROUTINE scm
2
3	USE flux_arp_mod_h
4	USE compbl_mod_h
5	USE clesphys_mod_h
6	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
7	USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
8	cwcon, clwcon, detr_therm, &
9	qsol, fevap, z0m, z0h, agesno, &
10	du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
11	falb_dir, falb_dif, &
12	ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
13	rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
14	solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, &
15	wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
16	wake_deltaq, wake_deltat, wake_s, awake_s, wake_dens, &
17	awake_dens, cv_gen, wake_cstar, &
18	zgam, zmax0, zmea, zpic, zsig, &
19	zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
20	ql_ancien, qs_ancien, qbs_ancien, &
21	cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
22	prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
23	u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
24
25
26	USE dimphy
27	USE surface_data, only : type_ocean,ok_veget
28	USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
29	pbl_surface_final
30	USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
31
32	USE infotrac ! new
33	USE control_mod
34	USE indice_sol_mod
35	USE phyaqua_mod
36	USE mod_1D_cases_read_std
37	USE print_control_mod, ONLY: lunout, prt_level
38	USE iniphysiq_mod, ONLY: iniphysiq
39	USE mod_const_mpi, ONLY: comm_lmdz
40	USE physiq_mod, ONLY: physiq
41	USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
42	preff, aps, bps, pseudoalt, scaleheight
43	USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
44	itau_dyn, itau_phy, start_time, year_len
45	USE phys_cal_mod, ONLY : year_len_phys_cal_mod => year_len
46	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_OUTPUTPHYSSCM
47
48	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
49	USE dimsoil_mod_h, ONLY: nsoilmx
50	USE yomcst_mod_h
51	USE tsoilnudge_mod_h
52	USE fcg_gcssold_mod_h
53	USE compar1d_mod_h
54	USE date_cas_mod_h
55	implicit none
56
57	!=====================================================================
58	! DECLARATIONS
59	!=====================================================================
60
61	!---------------------------------------------------------------------
62	! Externals
63	!---------------------------------------------------------------------
64	external fq_sat
65	real fq_sat
66
67	!---------------------------------------------------------------------
68	! Arguments d' initialisations de la physique (USER DEFINE)
69	!---------------------------------------------------------------------
70
71	integer, parameter :: ngrid=1
72	real :: zcufi = 1.
73	real :: zcvfi = 1.
74	real :: fnday
75	real :: day, daytime
76	real :: day1
77	real :: heure
78	integer :: jour
79	integer :: mois
80	integer :: an
81
82	!---------------------------------------------------------------------
83	! Declarations related to forcing and initial profiles
84	!---------------------------------------------------------------------
85
86	integer :: kmax = llm
87	integer llm700,nq1,nq2,iflag_1d_vert_adv
88	INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
89	real timestep, frac
90	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
91	real uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
92	real ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
93	real dqtdxls(nlev_max),dqtdyls(nlev_max)
94	real dqtdtls(nlev_max),thlpcar(nlev_max)
95	real qprof(nlev_max,nqmx)
96
97	! integer :: forcing_type
98	logical :: forcing_les = .false.
99	logical :: forcing_armcu = .false.
100	logical :: forcing_rico = .false.
101	logical :: forcing_radconv = .false.
102	logical :: forcing_toga = .false.
103	logical :: forcing_twpice = .false.
104	logical :: forcing_amma = .false.
105	logical :: forcing_dice = .false.
106	logical :: forcing_gabls4 = .false.
107
108	logical :: forcing_GCM2SCM = .false.
109	logical :: forcing_GCSSold = .false.
110	logical :: forcing_sandu = .false.
111	logical :: forcing_astex = .false.
112	logical :: forcing_fire = .false.
113	logical :: forcing_case = .false.
114	logical :: forcing_case2 = .false.
115	logical :: forcing_SCM = .false.
116
117	!flag forcings
118	logical :: nudge_wind=.true.
119	logical :: nudge_thermo=.false.
120	logical :: cptadvw=.true.
121
122
123	!=====================================================================
124	! DECLARATIONS FOR EACH CASE
125	!=====================================================================
126	!
127	INCLUDE "1D_decl_cases.h"
128	!
129	!---------------------------------------------------------------------
130	! Declarations related to nudging
131	!---------------------------------------------------------------------
132	integer :: nudge_max
133	parameter (nudge_max=9)
134	integer :: inudge_RHT=1
135	integer :: inudge_UV=2
136	logical :: nudge(nudge_max)
137	real :: t_targ(llm)
138	real :: rh_targ(llm)
139	real :: u_targ(llm)
140	real :: v_targ(llm)
141	!
142	!---------------------------------------------------------------------
143	! Declarations related to vertical discretization:
144	!---------------------------------------------------------------------
145	real :: pzero=1.e5
146	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
147	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
148
149	!---------------------------------------------------------------------
150	! Declarations related to variables
151	!---------------------------------------------------------------------
152
153	real :: phi(llm)
154	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
155	REAL rot(1, llm) ! relative vorticity, in s-1
156	real :: rlat_rad(1),rlon_rad(1)
157	real :: omega(llm),omega2(llm),rho(llm+1)
158	real :: ug(llm),vg(llm),fcoriolis
159	real :: sfdt, cfdt
160	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
161	real :: w_adv(llm),z_adv(llm)
162	real :: d_t_vert_adv(llm),d_u_vert_adv(llm),d_v_vert_adv(llm)
163	real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
164	real :: d_u_nudge(llm),d_v_nudge(llm)
165	real :: d_u_age(llm),d_v_age(llm)
166	real :: alpha
167	real :: ttt
168
169	REAL, ALLOCATABLE, DIMENSION(:,:):: q
170	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
171	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_vert_adv
172	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
173	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
174
175	!---------------------------------------------------------------------
176	! Initialization of surface variables
177	!---------------------------------------------------------------------
178	real :: run_off_lic_0(1)
179	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
180	real :: tsoil(1,nsoilmx,nbsrf)
181
182	!---------------------------------------------------------------------
183	! Call to phyredem
184	!---------------------------------------------------------------------
185	logical :: ok_writedem =.true.
186	real :: sollw_in = 0.
187	real :: solsw_in = 0.
188
189	!---------------------------------------------------------------------
190	! Call to physiq
191	!---------------------------------------------------------------------
192	logical :: firstcall=.true.
193	logical :: lastcall=.false.
194	real :: phis(1) = 0.0
195	real :: dpsrf(1)
196
197	!---------------------------------------------------------------------
198	! Initializations of boundary conditions
199	!---------------------------------------------------------------------
200	real, allocatable :: phy_nat (:) ! 0=ocean libre,1=land,2=glacier,3=banquise
201	real, allocatable :: phy_alb (:) ! Albedo land only (old value condsurf_jyg=0.3)
202	real, allocatable :: phy_sst (:) ! SST (will not be used; cf read_tsurf1d.F)
203	real, allocatable :: phy_bil (:) ! Ne sert que pour les slab_ocean
204	real, allocatable :: phy_rug (:) ! Longueur rugosite utilisee sur land only
205	real, allocatable :: phy_ice (:) ! Fraction de glace
206	real, allocatable :: phy_fter(:) ! Fraction de terre
207	real, allocatable :: phy_foce(:) ! Fraction de ocean
208	real, allocatable :: phy_fsic(:) ! Fraction de glace
209	real, allocatable :: phy_flic(:) ! Fraction de glace
210
211	!---------------------------------------------------------------------
212	! Fichiers et d'autres variables
213	!---------------------------------------------------------------------
214	integer :: k,l,i,it=1,mxcalc
215	integer :: nsrf
216	integer jcode
217	INTEGER read_climoz
218	integer :: it_end ! iteration number of the last call
219	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
220	data ecrit_slab_oc/-1/
221
222	! if flag_inhib_forcing = 0, tendencies of forcing are added
223	! <> 0, tendencies of forcing are not added
224	INTEGER :: flag_inhib_forcing = 0
225
226
227	print*,'VOUS ENTREZ DANS LE 1D FORMAT STANDARD'
228
229	!=====================================================================
230	! INITIALIZATIONS
231	!=====================================================================
232	du_phys(:)=0.
233	dv_phys(:)=0.
234	dt_phys(:)=0.
235	d_t_vert_adv(:)=0.
236	d_u_vert_adv(:)=0.
237	d_v_vert_adv(:)=0.
238	dt_cooling(:)=0.
239	d_t_adv(:)=0.
240	d_t_nudge(:)=0.
241	d_u_nudge(:)=0.
242	d_v_nudge(:)=0.
243	d_u_adv(:)=0.
244	d_v_adv(:)=0.
245	d_u_age(:)=0.
246	d_v_age(:)=0.
247
248
249	! Initialization of Common turb_forcing
250	dtime_frcg = 0.
251	Turb_fcg_gcssold=.false.
252	hthturb_gcssold = 0.
253	hqturb_gcssold = 0.
254
255
256
257
258	!---------------------------------------------------------------------
259	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
260	!---------------------------------------------------------------------
261	call conf_unicol
262	!Al1 moves this gcssold var from common fcg_gcssold to
263	Turb_fcg_gcssold = xTurb_fcg_gcssold
264	! --------------------------------------------------------------------
265	close(1)
266	write(,) 'lmdz1d.def lu => unicol.def'
267
268	forcing_SCM = .true.
269	year_ini_cas=1997
270	! It is possible that those parameters are run twice.
271	! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
272
273
274	call getin('anneeref',year_ini_cas)
275	call getin('dayref',day_deb)
276	mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
277	call getin('time_ini',heure_ini_cas)
278	print*,'NATURE DE LA SURFACE ',nat_surf
279
280	! Initialization of the logical switch for nudging
281
282	jcode = iflag_nudge
283	do i = 1,nudge_max
284	nudge(i) = mod(jcode,10) .ge. 1
285	jcode = jcode/10
286	enddo
287	! numerical scheme for vertical advection
288
289	iflag_1d_vert_adv=1
290	call getin('iflag_1d_vert_adv',iflag_1d_vert_adv)
291
292	!-----------------------------------------------------------------------
293	! Definition of the run
294	!-----------------------------------------------------------------------
295
296	call conf_gcm( 99, .TRUE. )
297
298	!-----------------------------------------------------------------------
299	allocate( phy_nat (year_len)) ! 0=ocean libre,1=land,2=glacier,3=banquise
300	phy_nat(:)=0.0
301	allocate( phy_alb (year_len)) ! Albedo land only (old value condsurf_jyg=0.3)
302	allocate( phy_sst (year_len)) ! SST (will not be used; cf read_tsurf1d.F)
303	allocate( phy_bil (year_len)) ! Ne sert que pour les slab_ocean
304	phy_bil(:)=1.0
305	allocate( phy_rug (year_len)) ! Longueur rugosite utilisee sur land only
306	allocate( phy_ice (year_len)) ! Fraction de glace
307	phy_ice(:)=0.0
308	allocate( phy_fter(year_len)) ! Fraction de terre
309	phy_fter(:)=0.0
310	allocate( phy_foce(year_len)) ! Fraction de ocean
311	phy_foce(:)=0.0
312	allocate( phy_fsic(year_len)) ! Fraction de glace
313	phy_fsic(:)=0.0
314	allocate( phy_flic(year_len)) ! Fraction de glace
315	phy_flic(:)=0.0
316
317
318	!-----------------------------------------------------------------------
319	! Choix du calendrier
320	! -------------------
321
322	! calend = 'earth_365d'
323	if (calend == 'earth_360d') then
324	call ioconf_calendar('360_day')
325	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
326	else if (calend == 'earth_365d') then
327	call ioconf_calendar('noleap')
328	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
329	else if (calend == 'earth_366d') then
330	call ioconf_calendar('all_leap')
331	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
332	else if (calend == 'gregorian') then
333	stop 'gregorian calend should not be used by normal user'
334	call ioconf_calendar('gregorian') ! not to be used by normal users
335	write(,)'CALENDRIER CHOISI: Gregorien'
336	else
337	write (,) 'ERROR : unknown calendar ', calend
338	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
339	endif
340	!-----------------------------------------------------------------------
341	!
342	!c Date :
343	! La date est supposee donnee sous la forme [annee, numero du jour dans
344	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
345	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
346	! Le numero du jour est dans "day". L heure est traitee separement.
347	! La date complete est dans "daytime" (l'unite est le jour).
348
349
350	if (nday>0) then
351	fnday=nday
352	else
353	fnday=-nday/float(day_step)
354	endif
355	print *,'fnday=',fnday
356	! start_time doit etre en FRACTION DE JOUR
357	start_time=time_ini/24.
358
359	annee_ref = anneeref
360	mois = 1
361	day_ref = dayref
362	heure = 0.
363	itau_dyn = 0
364	itau_phy = 0
365	call ymds2ju(annee_ref,mois,day_ref,heure,day)
366	day_ini = int(day)
367	day_end = day_ini + int(fnday)
368
369	! Convert the initial date to Julian day
370	day_ini_cas=day_deb
371	print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas
372	call ymds2ju &
373	& (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600 &
374	& ,day_ju_ini_cas)
375	print*,'time case 2',day_ini_cas,day_ju_ini_cas
376	daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation
377
378	! Print out the actual date of the beginning of the simulation :
379	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
380	print *,' Time of beginning : ', &
381	& year_print, month_print, day_print, sec_print
382
383	!---------------------------------------------------------------------
384	! Initialization of dimensions, geometry and initial state
385	!---------------------------------------------------------------------
386	! call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
387	! but we still need to initialize dimphy module (klon,klev,etc.) here.
388	call init_dimphy1D(1,llm)
389	call suphel
390	call init_infotrac
391
392	if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
393	allocate(q(llm,nqtot)) ; q(:,:)=0.
394	allocate(dq(llm,nqtot))
395	allocate(d_q_vert_adv(llm,nqtot))
396	allocate(d_q_adv(llm,nqtot))
397	allocate(d_q_nudge(llm,nqtot))
398
399	q(:,:) = 0.
400	dq(:,:) = 0.
401	d_q_vert_adv(:,:) = 0.
402	d_q_adv(:,:) = 0.
403	d_q_nudge(:,:) = 0.
404
405	! No ozone climatology need be read in this pre-initialization
406	! (phys_state_var_init is called again in physiq)
407	read_climoz = 0
408	nsw=6
409
410	call phys_state_var_init(read_climoz)
411
412	if (ngrid.ne.klon) then
413	print*,'stop in inifis'
414	print*,'Probleme de dimensions :'
415	print*,'ngrid = ',ngrid
416	print*,'klon = ',klon
417	stop
418	endif
419	!!!=====================================================================
420	!!! Feedback forcing values for Gateaux differentiation (al1)
421	!!!=====================================================================
422
423	qsol = qsolinp
424	qsurf = fq_sat(tsurf,psurf/100.)
425	beta_aridity(:,:) = beta_surf
426	day1= day_ini
427	time=daytime-day
428	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
429	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
430
431	!! mpl et jyg le 22/08/2012 :
432	!! pour que les cas a flux de surface imposes marchent
433	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
434	fsens=-wtsurfrcpdrho(1)
435	flat=-wqsurfrlvttrho(1)
436	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
437	ENDIF
438	print*,'Flux sol ',fsens,flat
439
440	! Vertical discretization and pressure levels at half and mid levels:
441
442	pa = 5e4
443	!! preff= 1.01325e5
444	preff = psurf
445	IF (ok_old_disvert) THEN
446	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
447	print *,'On utilise disvert0'
448	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
449	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
450	scaleheight=8.
451	pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
452	ELSE
453	call disvert()
454	print *,'On utilise disvert'
455	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
456	! Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
457	ENDIF
458
459	sig_s=presnivs/preff
460	plev =ap+bp*psurf
461	play = 0.5*(plev(1:llm)+plev(2:llm+1))
462	zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles.
463
464	IF (forcing_type .eq. 59) THEN
465	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
466	write(,) '***********************'
467	do l = 1, llm
468	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
469	if (trouve_700 .and. play(l).le.70000) then
470	llm700=l
471	print *,'llm700,play=',llm700,play(l)/100.
472	trouve_700= .false.
473	endif
474	enddo
475	write(,) '***********************'
476	ENDIF
477
478	!
479	!=====================================================================
480	! EVENTUALLY, READ FORCING DATA :
481	!=====================================================================
482
483	INCLUDE "1D_read_forc_cases.h"
484
485	if (forcing_GCM2SCM) then
486	write (,) 'forcing_GCM2SCM not yet implemented'
487	stop 'in initialization'
488	endif ! forcing_GCM2SCM
489
490
491	!=====================================================================
492	! Initialisation de la physique :
493	!=====================================================================
494
495	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
496	!
497	! day_step, iphysiq lus dans gcm.def ci-dessus
498	! timestep: calcule ci-dessous from rday et day_step
499	! ngrid=1
500	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
501	! rday: defini dans suphel.F (86400.)
502	! day_ini: lu dans run.def (dayref)
503	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
504	! airefi,zcufi,zcvfi initialises au debut de ce programme
505	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
506
507
508	day_step = float(nsplit_phys)*day_step/float(iphysiq)
509	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
510	timestep =rday/day_step
511	dtime_frcg = timestep
512	!
513	zcufi=airefi
514	zcvfi=airefi
515	!
516	rlat_rad(1)=xlat*rpi/180.
517	rlon_rad(1)=xlon*rpi/180.
518
519	! iniphysiq will call iniaqua who needs year_len from phys_cal_mod
520	year_len_phys_cal_mod=year_len
521
522	! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
523	! e.g. for cell boundaries, which are meaningless in 1D; so pad these
524	! with '0.' when necessary
525
526	call iniphysiq(iim,jjm,llm, &
527	1,comm_lmdz, &
528	rday,day_ini,timestep, &
529	(/rlat_rad(1),0./),(/0./), &
530	(/0.,0./),(/rlon_rad(1),0./), &
531	(/ (/airefi,0./),(/0.,0./) /), &
532	(/zcufi,0.,0.,0./), &
533	(/zcvfi,0./), &
534	ra,rg,rd,rcpd,1)
535	print*,'apres iniphysiq'
536
537	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
538	co2_ppm= 330.0
539	solaire=1370.0
540
541	! Ecriture du startphy avant le premier appel a la physique.
542	! On le met juste avant pour avoir acces a tous les champs
543
544	if (ok_writedem) then
545
546	!--------------------------------------------------------------------------
547	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
548	! need : qsol fder snow qsurf evap rugos agesno ftsoil
549	!--------------------------------------------------------------------------
550
551	type_ocean = "force"
552	run_off_lic_0(1) = restart_runoff
553	call fonte_neige_init(run_off_lic_0)
554
555	fder=0.
556	snsrf(1,:)=snowmass ! masse de neige des sous surface
557	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
558	fevap=0.
559	z0m(1,:)=rugos ! couverture de neige des sous surface
560	z0h(1,:)=rugosh ! couverture de neige des sous surface
561	agesno = xagesno
562	tsoil(:,:,:)=tsurf
563	!-----------------------------------------------------------------------
564	call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
565
566	!------------------ prepare limit conditions for limit.nc -----------------
567	!-- Ocean force
568
569	print*,'avant phyredem'
570	pctsrf(1,:)=0.
571	if (nat_surf.eq.0.) then
572	pctsrf(1,is_oce)=1.
573	pctsrf(1,is_ter)=0.
574	pctsrf(1,is_lic)=0.
575	pctsrf(1,is_sic)=0.
576	else if (nat_surf .eq. 1) then
577	pctsrf(1,is_oce)=0.
578	pctsrf(1,is_ter)=1.
579	pctsrf(1,is_lic)=0.
580	pctsrf(1,is_sic)=0.
581	else if (nat_surf .eq. 2) then
582	pctsrf(1,is_oce)=0.
583	pctsrf(1,is_ter)=0.
584	pctsrf(1,is_lic)=1.
585	pctsrf(1,is_sic)=0.
586	else if (nat_surf .eq. 3) then
587	pctsrf(1,is_oce)=0.
588	pctsrf(1,is_ter)=0.
589	pctsrf(1,is_lic)=0.
590	pctsrf(1,is_sic)=1.
591
592	end if
593
594
595	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf &
596	& ,pctsrf(1,is_oce),pctsrf(1,is_ter)
597
598	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
599	zpic = zpicinp
600	ftsol=tsurf
601	falb_dir=albedo
602	falb_dif=albedo
603	rugoro=rugos
604	t_ancien(1,:)=temp(:)
605	q_ancien(1,:)=q(:,1)
606	ql_ancien = 0.
607	qs_ancien = 0.
608	prlw_ancien = 0.
609	prsw_ancien = 0.
610	prw_ancien = 0.
611	IF ( ok_bs ) THEN
612	qbs_ancien = 0.
613	prbsw_ancien = 0.
614	ENDIF
615	IF ( ok_ice_supersat ) THEN
616	cf_ancien = 0.
617	qvc_ancien = 0.
618	cwcon = 0.
619	ENDIF
620	IF ( ok_plane_contrail ) THEN
621	cfa_ancien = 0.
622	pcf_ancien = 0.
623	qva_ancien = 0.
624	qia_ancien = 0.
625	ENDIF
626	rain_fall=0.
627	snow_fall=0.
628	solsw=0.
629	solswfdiff=0.
630	sollw=0.
631	sollwdown=rsigmatsurf*4
632	radsol=0.
633	rnebcon=0.
634	ratqs=0.
635	clwcon=0.
636	zmax0 = 0.
637	zmea=zsurf
638	zstd=0.
639	zsig=0.
640	zgam=0.
641	zval=0.
642	zthe=0.
643	sig1=0.
644	w01=0.
645	!
646	wake_deltaq = 0.
647	wake_deltat = 0.
648	wake_delta_pbl_TKE(:,:,:) = 0.
649	delta_tsurf = 0.
650	wake_fip = 0.
651	wake_pe = 0.
652	wake_s = 0.
653	awake_s = 0.
654	wake_dens = 0.
655	awake_dens = 0.
656	cv_gen = 0.
657	wake_cstar = 0.
658	ale_bl = 0.
659	ale_bl_trig = 0.
660	alp_bl = 0.
661	IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
662	IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
663	entr_therm = 0.
664	detr_therm = 0.
665	f0 = 0.
666	fm_therm = 0.
667	u_ancien(1,:)=u(:)
668	v_ancien(1,:)=v(:)
669
670	u10m=0.
671	v10m=0.
672	ale_wake=0.
673	ale_bl_stat=0.
674	ratqs_inter_(:,:)= 0.002
675
676	!------------------------------------------------------------------------
677	! Make file containing restart for the physics (startphy.nc)
678	!
679	! NB: List of the variables to be written by phyredem (via put_field):
680	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
681	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
682	! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
683	! radsol,solsw,solswfdiff,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
684	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
685	! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
686	! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
687	! wake_deltat,wake_deltaq,wake_s,awake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
688	! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
689	!
690	! NB2: The content of the startphy.nc file depends on some flags defined in
691	! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
692	! to be set at some arbitratry convenient values.
693	!------------------------------------------------------------------------
694	!Al1 =============== restart option ======================================
695	iflag_physiq=0
696	call getin('iflag_physiq',iflag_physiq)
697
698	if (.not.restart) then
699	iflag_pbl = 5
700	call phyredem ("startphy.nc")
701	else
702	! (desallocations)
703	print*,'callin surf final'
704	call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
705	print*,'after surf final'
706	CALL fonte_neige_final(run_off_lic_0)
707	endif
708
709	ok_writedem=.false.
710	print*,'apres phyredem'
711
712	endif ! ok_writedem
713
714	!------------------------------------------------------------------------
715	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
716	! --------------------------------------------------
717	! NB: List of the variables to be written in limit.nc
718	! (by writelim.F, subroutine of 1DUTILS.h):
719	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
720	! phy_fter,phy_foce,phy_flic,phy_fsic)
721	!------------------------------------------------------------------------
722	do i=1,year_len
723	phy_nat(i) = nat_surf
724	phy_alb(i) = albedo
725	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
726	phy_rug(i) = rugos
727	phy_fter(i) = pctsrf(1,is_ter)
728	phy_foce(i) = pctsrf(1,is_oce)
729	phy_fsic(i) = pctsrf(1,is_sic)
730	phy_flic(i) = pctsrf(1,is_lic)
731	enddo
732
733	! fabrication de limit.nc
734	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug, &
735	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
736
737
738	call phys_state_var_end
739	if (restart) then
740	print*,'call to restart dyn 1d'
741	Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs, &
742	& u,v,temp,q,omega2)
743
744	print*,'fnday,annee_ref,day_ref,day_ini', &
745	& fnday,annee_ref,day_ref,day_ini
746	day = day_ini
747	day_end = day_ini + nday
748	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
749
750	! Print out the actual date of the beginning of the simulation :
751	call ju2ymds(daytime, an, mois, jour, heure)
752	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
753
754	day = int(daytime)
755	time=daytime-day
756
757	print,'*** intialised fields from restart1dyn *****'
758	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
759	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
760	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
761	! raz for safety
762	do l=1,llm
763	d_q_vert_adv(l,:) = 0.
764	enddo
765	endif
766	!====================== end restart =================================
767	IF (ecrit_slab_oc.eq.1) then
768	open(97,file='div_slab.dat',STATUS='UNKNOWN')
769	elseif (ecrit_slab_oc.eq.0) then
770	open(97,file='div_slab.dat',STATUS='OLD')
771	endif
772	!
773	!=====================================================================
774	IF (CPPKEY_OUTPUTPHYSSCM) THEN
775	CALL iophys_ini(timestep)
776	END IF
777
778	!=====================================================================
779	! START OF THE TEMPORAL LOOP :
780	!=====================================================================
781
782	it_end = nint(fnday*day_step)
783	do while(it.le.it_end)
784
785	if (prt_level.ge.1) then
786	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
787	& it,day,time,it_end,day_step
788	print*,'PAS DE TEMPS ',timestep
789	endif
790	if (it.eq.it_end) lastcall=.True.
791
792	!---------------------------------------------------------------------
793	! Interpolation of forcings in time and onto model levels
794	!---------------------------------------------------------------------
795
796	INCLUDE "1D_interp_cases.h"
797
798	!---------------------------------------------------------------------
799	! Geopotential :
800	!---------------------------------------------------------------------
801	phis(1)=zsurf*RG
802
803	! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
804	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
805
806	do l = 1, llm-1
807	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
808	& (play(l)-play(l+1))/(play(l)+play(l+1))
809	enddo
810
811	!---------------------------------------------------------------------
812	! Vertical advection
813	!---------------------------------------------------------------------
814
815	IF ( forc_w+forc_omega > 0 ) THEN
816
817	IF ( forc_w == 1 ) THEN
818	w_adv=w_mod_cas
819	z_adv=phi/RG
820	ELSE
821	w_adv=omega
822	z_adv=play
823	ENDIF
824
825	! calculation of potential temperature for the advection
826	teta=temp(pzero/play)*rkappa
827
828	! vertical tendencies computed as d X / d t = -W d X / d z
829
830	IF (iflag_1d_vert_adv .EQ. 0) THEN
831
832	! old centered numerical scheme
833	do l=2,llm-1
834	d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
835	d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
836	! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
837	d_t_vert_adv(l)=-w_adv(l)(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))*rkappa
838	d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1))
839	enddo
840
841
842	ELSE IF (iflag_1d_vert_adv .EQ. 1) THEN
843	! upstream numerical scheme
844
845	do l=2,llm-1
846	IF ( ( ( forc_w .EQ. 1 ) .AND. ( w_adv(l) .GT. 0. ) ) .OR. &
847	( ( forc_w .NE. 1 ) .AND. ( w_adv(l) .LT. 0. ) ) ) THEN
848	d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1))
849	d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1))
850	! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
851	d_t_vert_adv(l)=-w_adv(l)(teta(l)-teta(l-1))/(z_adv(l)-z_adv(l-1)) / (pzero/play(l))*rkappa
852	d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1))
853	ELSE
854	d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l))
855	d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l))
856	! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
857	d_t_vert_adv(l)=-w_adv(l)(teta(l+1)-teta(l))/(z_adv(l+1)-z_adv(l)) / (pzero/play(l))*rkappa
858	d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l))
859	ENDIF
860	enddo
861
862
863	ENDIF ! numerical scheme for vertical advection
864
865
866	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
867	u(1:mxcalc)=u(1:mxcalc) + timestep * d_u_vert_adv(1:mxcalc)
868	v(1:mxcalc)=v(1:mxcalc) + timestep * d_v_vert_adv(1:mxcalc)
869	q(1:mxcalc,:)=q(1:mxcalc,:) + timestep * d_q_vert_adv(1:mxcalc,:)
870	temp(1:mxcalc)=temp(1:mxcalc) + timestep * d_t_vert_adv(1:mxcalc)
871	teta=temp(pzero/play)*rkappa
872	ENDIF
873
874	ENDIF
875
876	!---------------------------------------------------------------------
877	! Listing output for debug prt_level>=1
878	!---------------------------------------------------------------------
879	if (prt_level>=1) then
880	print *,' avant physiq : -------- day time ',day,time
881	write(,) 'firstcall,lastcall,phis', &
882	& firstcall,lastcall,phis
883	end if
884	if (prt_level>=5) then
885	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l', &
886	& 'presniv','plev','play','phi'
887	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l, &
888	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
889	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l', &
890	& 'presniv','u','v','temp','q1','q2','omega2'
891	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l, &
892	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
893	endif
894
895	!---------------------------------------------------------------------
896	! Call physiq :
897	!---------------------------------------------------------------------
898	call physiq(ngrid,llm, &
899	firstcall,lastcall,timestep, &
900	plev,play,phi,phis,presnivs, &
901	u,v, rot, temp,q,omega2, &
902	du_phys,dv_phys,dt_phys,dq,dpsrf)
903	firstcall=.false.
904
905	!---------------------------------------------------------------------
906	! Listing output for debug
907	!---------------------------------------------------------------------
908	if (prt_level>=5) then
909	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l', &
910	& 'presniv','plev','play','phi'
911	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l, &
912	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
913	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l', &
914	& 'presniv','u','v','temp','q1','q2','omega2'
915	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l, &
916	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
917	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l', &
918	& 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
919	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l, &
920	& presnivs(l),86400du_phys(l),86400dv_phys(l), &
921	& 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
922	write(,) 'dpsrf',dpsrf
923	endif
924	!---------------------------------------------------------------------
925	! Add physical tendencies :
926	!---------------------------------------------------------------------
927
928	fcoriolis=2.sin(rpixlat/180.)*romega
929
930	IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
931	fcoriolis, xlat,mxcalc
932
933	!---------------------------------------------------------------------
934	! Geostrophic forcing
935	!---------------------------------------------------------------------
936
937	IF ( forc_geo == 0 ) THEN
938	d_u_age(1:mxcalc)=0.
939	d_v_age(1:mxcalc)=0.
940	ELSE
941	sfdt = sin(0.5fcoriolistimestep)
942	cfdt = cos(0.5fcoriolistimestep)
943
944	d_u_age(1:mxcalc)= -2.sfdt/timestep &
945	& (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - &
946	& cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
947	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
948	!
949	d_v_age(1:mxcalc)= -2.sfdt/timestep &
950	& (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + &
951	& sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
952	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
953	ENDIF
954	!
955	!---------------------------------------------------------------------
956	! Nudging
957	!---------------------------------------------------------------------
958	d_t_nudge(:) = 0.
959	d_u_nudge(:) = 0.
960	d_v_nudge(:) = 0.
961	d_q_nudge(:,:) = 0.
962
963	DO l=1,llm
964
965	IF (nudging_u .LT. 0) THEN
966
967	d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))*invtau_u_nudg_mod_cas(l)
968
969	ELSE
970
971	IF ( play(l) < p_nudging_u .AND. nint(nudging_u) /= 0 ) &
972	& d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))/nudging_u
973
974	ENDIF
975
976
977	IF (nudging_v .LT. 0) THEN
978
979	d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))*invtau_v_nudg_mod_cas(l)
980
981	ELSE
982
983
984	IF ( play(l) < p_nudging_v .AND. nint(nudging_v) /= 0 ) &
985	& d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))/nudging_v
986
987	ENDIF
988
989
990	IF (nudging_t .LT. 0) THEN
991
992	d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))*invtau_temp_nudg_mod_cas(l)
993
994	ELSE
995
996
997	IF ( play(l) < p_nudging_t .AND. nint(nudging_t) /= 0 ) &
998	& d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))/nudging_t
999
1000	ENDIF
1001
1002
1003	IF (nudging_qv .LT. 0) THEN
1004
1005	d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))*invtau_qv_nudg_mod_cas(l)
1006
1007	ELSE
1008
1009	IF ( play(l) < p_nudging_qv .AND. nint(nudging_qv) /= 0 ) &
1010	& d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv
1011
1012	ENDIF
1013
1014	ENDDO
1015
1016	!-----------------------------------------------------------
1017	! horizontal forcings (advection) and nudging
1018	!-----------------------------------------------------------
1019
1020	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
1021
1022	u(1:mxcalc)=u(1:mxcalc) + timestep*( &
1023	& du_phys(1:mxcalc) &
1024	& +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) &
1025	& +d_u_nudge(1:mxcalc) )
1026	v(1:mxcalc)=v(1:mxcalc) + timestep*( &
1027	& dv_phys(1:mxcalc) &
1028	& +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) &
1029	& +d_v_nudge(1:mxcalc) )
1030	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( &
1031	& dq(1:mxcalc,:) &
1032	& +d_q_adv(1:mxcalc,:) &
1033	& +d_q_nudge(1:mxcalc,:) )
1034
1035
1036	temp(1:mxcalc)=temp(1:mxcalc)+timestep*( &
1037	& dt_phys(1:mxcalc) &
1038	& +d_t_adv(1:mxcalc) &
1039	& +d_t_nudge(1:mxcalc) &
1040	& +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
1041
1042	!---------------------------------------------------------------------
1043	! Optional outputs
1044	!---------------------------------------------------------------------
1045
1046	IF (CPPKEY_OUTPUTPHYSSCM) THEN
1047	CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
1048	CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv)
1049	CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv)
1050	CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv)
1051	CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1))
1052	CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas)
1053	CALL iophys_ecrit('u',klev,'u','m/s',u)
1054	CALL iophys_ecrit('unud',klev,'unud','m/s',u_nudg_mod_cas)
1055	CALL iophys_ecrit('v',klev,'v','m/s',v)
1056	CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas)
1057	CALL iophys_ecrit('temp',klev,'temp','K',temp)
1058	CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas)
1059	CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge)
1060	CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',100086400d_q_nudge(:,1))
1061	END IF
1062
1063
1064
1065	! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
1066
1067	teta=temp(pzero/play)*rkappa
1068
1069	!---------------------------------------------------------------------
1070	! Nudge soil temperature if requested
1071	!---------------------------------------------------------------------
1072
1073	IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1074	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) &
1075	& -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
1076	ENDIF
1077
1078
1079	END IF ! end if in case tendency should be added
1080
1081	!---------------------------------------------------------------------
1082	! Air temperature :
1083	!---------------------------------------------------------------------
1084	if (lastcall) then
1085	print*,'Pas de temps final ',it
1086	call ju2ymds(daytime, an, mois, jour, heure)
1087	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1088	endif
1089
1090	! incremente day time
1091	daytime = daytime+1./day_step
1092	day = int(daytime+0.1/day_step)
1093	time = time_ini/24.+real(mod(it,day_step))/day_step
1094	it=it+1
1095
1096	enddo
1097
1098	if (ecrit_slab_oc.ne.-1) close(97)
1099
1100	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1101	! ---------------------------------------------------------------------------
1102	call dyn1dredem("restart1dyn.nc", &
1103	& plev,play,phi,phis,presnivs, &
1104	& u,v,temp,q,omega2)
1105
1106	CALL abort_gcm ('lmdz1d ','The End ',0)
1107
1108	END SUBROUTINE scm

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

Download in other formats:

Original Format