Context Navigation

scm.F90 @ 3541

Last change on this file since 3541 was 3541, checked in by fhourdin, 5 years ago

Gros nettoyage en cours sur le 1D.
Le nouveau lmdz1d.F90 est une coquille vide qui choisit entre
old_lmdz1d.F90 (l'ancien lmdz1d.F90) et scm.F90 (le nouveau au format standard).
Plusieur fichiers sont renommés old_truc, le truc étant au format standard,
nettoyé des anciens formats.
Le 1DUTILS.h est lui même séparé entre des routines génériques venant remplacer
notamment des routines de dyn3d (la vocation d'origine de 1DUTILS.h) et
les routiles de lecture spécifiques mises dans old_1DUTILS.h
On perdra un peu de l'utilister de svn au moment de cette grosse bascule.
Mais les old_ sont faits pour ne plus bouger, et les versions standard
sont en pleine évolution.
Fredho

File size: 40.9 KB

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

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

Download in other formats:

Original Format