Context Navigation

lmdz1d.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706

File size: 28.7 KB

Line
1	PROGRAM lmdz1d
2
3	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
4	use phys_state_var_mod
5	use comgeomphy
6	use dimphy
7	use surface_data, only : type_ocean,ok_veget
8	use pbl_surface_mod, only : pbl_surface_init, pbl_surface_final
9	use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
10
11	use infotrac ! new
12	use control_mod
13
14	implicit none
15	#include "dimensions.h"
16	#include "YOMCST.h"
17	#include "temps.h"
18	!!#include "control.h"
19	#include "iniprint.h"
20	#include "clesphys.h"
21	#include "dimsoil.h"
22	#include "indicesol.h"
23
24	#include "comvert.h"
25	#include "compar1d.h"
26	#include "flux_arp.h"
27	#include "fcg_gcssold.h"
28	!!!#include "fbforcing.h"
29
30	!=====================================================================
31	! DECLARATIONS
32	!=====================================================================
33
34	!---------------------------------------------------------------------
35	! Externals
36	!---------------------------------------------------------------------
37	external fq_sat
38	real fq_sat
39
40	!---------------------------------------------------------------------
41	! Arguments d' initialisations de la physique (USER DEFINE)
42	!---------------------------------------------------------------------
43
44	integer, parameter :: ngrid=1
45	real :: zcufi = 1.
46	real :: zcvfi = 1.
47
48	!- real :: nat_surf
49	!- logical :: ok_flux_surf
50	!- real :: fsens
51	!- real :: flat
52	!- real :: tsurf
53	!- real :: rugos
54	!- real :: qsol(1:2)
55	!- real :: qsurf
56	!- real :: psurf
57	!- real :: zsurf
58	!- real :: albedo
59	!-
60	!- real :: time = 0.
61	!- real :: time_ini
62	!- real :: xlat
63	!- real :: xlon
64	!- real :: wtsurf
65	!- real :: wqsurf
66	!- real :: restart_runoff
67	!- real :: xagesno
68	!- real :: qsolinp
69	!- real :: zpicinp
70	!-
71	real :: fnday
72	real :: day, daytime
73	real :: day1
74	real :: heure
75	integer :: jour
76	integer :: mois
77	integer :: an
78
79	!
80	real :: paire = 1. ! aire de la maille
81	!** common /flux_arp/fsens,flat,ok_flux_surf
82
83	!---------------------------------------------------------------------
84	! Declarations related to forcing and initial profiles
85	!---------------------------------------------------------------------
86
87	integer :: kmax = llm
88	integer nlev_max
89	parameter (nlev_max = 100)
90	real timestep, frac, timeit
91	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max),
92	. uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
93	. ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
94	. dqtdxls(nlev_max),dqtdyls(nlev_max),
95	. dqtdtls(nlev_max),thlpcar(nlev_max)
96
97	real :: fff
98	c integer :: forcing_type
99	logical :: forcing_les = .false.
100	logical :: forcing_armcu = .false.
101	logical :: forcing_rico = .false.
102	logical :: forcing_radconv = .false.
103	logical :: forcing_toga = .false.
104	logical :: forcing_twpice = .false.
105	logical :: forcing_GCM2SCM = .false.
106	logical :: forcing_GCSSold = .false.
107	integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
108	! (cf read_tsurf1d.F)
109
110	!vertical advection computation
111	real d_t_z(llm), d_q_z(llm)
112	real d_t_dyn_z(llm), d_q_dyn_z(llm)
113	real zz(llm)
114	real zfact
115
116	!flag forcings
117	logical :: nudge_wind=.true.
118	logical :: nudge_thermo=.false.
119	logical :: cptadvw=.true.
120	!=====================================================================
121	! DECLARATIONS FOR EACH CASE
122	!=====================================================================
123	!
124	#include "1D_decl_cases.h"
125	!
126	!---------------------------------------------------------------------
127	! Declarations related to vertical discretization:
128	!---------------------------------------------------------------------
129	real :: pzero=1.e5
130	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
131	real :: playd(llm),zlayd(llm)
132
133	!---------------------------------------------------------------------
134	! Declarations related to variables
135	!---------------------------------------------------------------------
136
137	integer :: iq
138	real :: phi(llm)
139	real :: teta(llm),temp(llm),u(llm),v(llm)
140	real :: omega(llm+1),omega2(llm),rho(llm+1)
141	real :: ug(llm),vg(llm),fcoriolis
142	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
143	real :: du_dyn(llm),dv_dyn(llm),dt_dyn(llm)
144	real :: dt_cooling(llm),d_t_cool(llm),d_th_adv(llm)
145	real :: dq_cooling(llm),d_q_cool(llm)
146	real :: tmpvar(llm)
147	real :: alpha
148
149	REAL, ALLOCATABLE, DIMENSION(:,:):: q
150	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
151	REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
152	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
153
154	!---------------------------------------------------------------------
155	! Initialization of surface variables
156	!---------------------------------------------------------------------
157	real :: run_off_lic_0(1)
158	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
159	real :: evap(1,nbsrf),frugs(1,nbsrf)
160	real :: tsoil(1,nsoilmx,nbsrf)
161	real :: agesno(1,nbsrf)
162
163	!---------------------------------------------------------------------
164	! Call to phyredem
165	!---------------------------------------------------------------------
166	logical :: ok_writedem =.true.
167
168	!---------------------------------------------------------------------
169	! Call to physiq
170	!---------------------------------------------------------------------
171	integer, parameter :: longcles=20
172	logical :: firstcall=.true.
173	logical :: lastcall=.false.
174	real :: phis = 0.0
175	real :: clesphy0(longcles) = 0.0
176	real :: dpsrf
177
178	!---------------------------------------------------------------------
179	! Initializations of boundary conditions
180	!---------------------------------------------------------------------
181	integer, parameter :: yd = 360
182	real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
183	real :: phy_alb (yd) ! Albedo land only (old value condsurf_jyg=0.3)
184	real :: phy_sst (yd) ! SST (will not be used; cf read_tsurf1d.F)
185	real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
186	real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
187	real :: phy_ice (yd) = 0.0 ! Fraction de glace
188	real :: phy_fter(yd) = 0.0 ! Fraction de terre
189	real :: phy_foce(yd) = 0.0 ! Fraction de ocean
190	real :: phy_fsic(yd) = 0.0 ! Fraction de glace
191	real :: phy_flic(yd) = 0.0 ! Fraction de glace
192
193	!---------------------------------------------------------------------
194	! Fichiers et d'autres variables
195	!---------------------------------------------------------------------
196	real ttt
197	integer :: ierr,k,l,i,it=1,mxcalc
198	integer jjmp1
199	parameter (jjmp1=jjm+1-1/jjm)
200	INTEGER nbteta
201	PARAMETER(nbteta=3)
202	REAL dudyn(iim+1,jjmp1,llm)
203	REAL PVteta(1,nbteta)
204	INTEGER read_climoz
205	!Al1
206	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
207	data ecrit_slab_oc/-1/
208
209	!=====================================================================
210	! INITIALIZATIONS
211	!=====================================================================
212	! Initialization of Common turb_forcing
213	dtime_frcg = 0.
214	Turb_fcg_gcssold=.false.
215	hthturb_gcssold = 0.
216	hqturb_gcssold = 0.
217
218	!---------------------------------------------------------------------
219	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
220	!---------------------------------------------------------------------
221	cAl1
222	call conf_unicol(99)
223	cAl1 moves this gcssold var from common fcg_gcssold to
224	Turb_fcg_gcssold = xTurb_fcg_gcssold
225	c --------------------------------------------------------------------
226	close(1)
227	cAl1
228	write(,) 'lmdz1d.def lu => unicol.def'
229
230	! forcing_type defines the way the SCM is forced:
231	!forcing_type = 0 ==> forcing_les = .true.
232	! initial profiles from file prof.inp.001
233	! no forcing by LS convergence ;
234	! surface temperature imposed ;
235	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
236	!forcing_type = 1 ==> forcing_radconv = .true.
237	! idem forcing_type = 0, but the imposed radiative cooling
238	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
239	! then there is no radiative cooling at all)
240	!forcing_type = 2 ==> forcing_toga = .true.
241	! initial profiles from TOGA-COARE IFA files
242	! LS convergence and SST imposed from TOGA-COARE IFA files
243	!forcing_type = 3 ==> forcing_GCM2SCM = .true.
244	! initial profiles from the GCM output
245	! LS convergence imposed from the GCM output
246	!forcing_type = 4 ==> forcing_twpice = .true.
247	! initial profiles from TWP-ICE cdf file
248	! LS convergence, omega and SST imposed from TWP-ICE files
249	!forcing_type = 5 ==> forcing_rico = .true.
250	! initial profiles from RICO files
251	! LS convergence imposed from RICO files
252	!forcing_type = 40 ==> forcing_GCSSold = .true.
253	! initial profile from GCSS file
254	! LS convergence imposed from GCSS file
255	!forcing_type = 61 ==> forcing_armcu = .true.
256	! initial profile from arm_cu file
257	! LS convergence imposed from arm_cu file
258	!
259	if (forcing_type .eq.0) THEN
260	forcing_les = .true.
261	elseif (forcing_type .eq.1) THEN
262	forcing_radconv = .true.
263	elseif (forcing_type .eq.2) THEN
264	forcing_toga = .true.
265	elseif (forcing_type .eq.3) THEN
266	forcing_GCM2SCM = .true.
267	elseif (forcing_type .eq.4) THEN
268	forcing_twpice = .true.
269	elseif (forcing_type .eq.5) THEN
270	forcing_rico = .true.
271	elseif (forcing_type .eq.40) THEN
272	forcing_GCSSold = .true.
273	elseif (forcing_type .eq.61) THEN
274	forcing_armcu = .true.
275	IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
276	else
277	write (,) 'ERROR : unknown forcing_type ', forcing_type
278	stop 'Forcing_type should be 0,1,2,3 or 40'
279	ENDIF
280	print*,"forcing type=",forcing_type
281
282	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
283	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
284	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
285	! through the common sst_forcing.
286
287	type_ts_forcing = 0
288	if (forcing_toga) type_ts_forcing = 1
289
290	!---------------------------------------------------------------------
291	! Definition of the run
292	!---------------------------------------------------------------------
293
294	call conf_gcm( 99, .TRUE. , clesphy0 )
295	c-----------------------------------------------------------------------
296	c Choix du calendrier
297	c -------------------
298
299	c calend = 'earth_365d'
300	if (calend == 'earth_360d') then
301	call ioconf_calendar('360d')
302	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
303	else if (calend == 'earth_365d') then
304	call ioconf_calendar('noleap')
305	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
306	else if (calend == 'earth_366d') then
307	call ioconf_calendar('all_leap')
308	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
309	else if (calend == 'gregorian') then
310	call ioconf_calendar('gregorian') ! not to be used by normal users
311	write(,)'CALENDRIER CHOISI: Gregorien'
312	else
313	write (,) 'ERROR : unknown calendar ', calend
314	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
315	endif
316	c-----------------------------------------------------------------------
317	c
318	cc Date :
319	c La date est supposee donnee sous la forme [annee, numero du jour dans
320	c l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
321	c On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
322	c Le numero du jour est dans "day". L heure est traitee separement.
323	c La date complete est dans "daytime" (l'unite est le jour).
324	fnday=nday
325	c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
326	IF(forcing_type .EQ. 61) fnday=53100./86400.
327	annee_ref = anneeref
328	mois = 1
329	day_ref = dayref
330	heure = 0.
331	itau_dyn = 0
332	itau_phy = 0
333	call ymds2ju(annee_ref,mois,day_ref,heure,day)
334	day_ini = day
335	day_end = day_ini + nday
336	! Convert the initial date of Toga-Coare to Julian day
337	call ymds2ju
338	$ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
339
340	! Convert the initial date of TWPICE to Julian day
341	call ymds2ju
342	$ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
343	$ ,day_ju_ini_twpi)
344
345	! Convert the initial date of Arm_cu to Julian day
346	call ymds2ju
347	$ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
348	$ ,day_ju_ini_armcu)
349
350	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
351	! Print out the actual date of the beginning of the simulation :
352	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
353	print *,' Time of beginning : ',
354	$ year_print, month_print, day_print, sec_print
355
356	!---------------------------------------------------------------------
357	! Initialization of dimensions, geometry and initial state
358	!---------------------------------------------------------------------
359	call init_phys_lmdz(1,1,llm,1,(/1/))
360	call suphel
361	call initcomgeomphy
362	call infotrac_init
363
364	allocate(q(llm,nqtot))
365	allocate(dq(llm,nqtot))
366	allocate(dq_dyn(llm,nqtot))
367	allocate(d_q_adv(llm,nqtot))
368
369	c
370	c No ozone climatology need be read in this pre-initialization
371	c (phys_state_var_init is called again in physiq)
372	read_climoz = 0
373	c
374	call phys_state_var_init(read_climoz)
375
376	if (ngrid.ne.klon) then
377	print*,'stop in inifis'
378	print*,'Probleme de dimensions :'
379	print*,'ngrid = ',ngrid
380	print*,'klon = ',klon
381	stop
382	endif
383	!!!=====================================================================
384	!!! Feedback forcing values for Gateaux differentiation (al1)
385	!!!=====================================================================
386	!!! Surface Planck forcing bracketing call radiation
387	!! surf_Planck = 0.
388	!! surf_Conv = 0.
389	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
390	!!! a mettre dans le lmdz1d.def ou autre
391	!!
392	!!
393	qsol = qsolinp
394	qsurf = fq_sat(tsurf,psurf/100.)
395	rlat=xlat
396	rlon=xlon
397	day1= day_ini
398	time=daytime-day
399	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
400	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
401
402	!
403	!! mpl et jyg le 22/08/2012 :
404	!! pour que les cas a flux de surface imposes marchent
405	IF(.NOT.ok_flux_surf) THEN
406	fsens=-wtsurfrcpdrho(1)
407	flat=-wqsurfrlvttrho(1)
408	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
409	ENDIF
410	!! ok_flux_surf=.false.
411	!! fsens=-wtsurfrcpdrho(1)
412	!! flat=-wqsurfrlvttrho(1)
413	!!!!
414
415	! Vertical discretization and pressure levels at half and mid levels:
416
417	pa = 5e4
418	!! preff= 1.01325e5
419	preff = psurf
420	call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
421	sig_s=presnivs/preff
422	plev =ap+bp*psurf
423	play = 0.5*(plev(1:llm)+plev(2:llm+1))
424	ccc zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles
425
426	write(,) '***********************'
427	do l = 1, llm
428	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
429	enddo
430	write(,) '***********************'
431
432	c
433	!=====================================================================
434	! EVENTUALLY, READ FORCING DATA :
435	!=====================================================================
436
437	#include "1D_read_forc_cases.h"
438
439	if (forcing_GCM2SCM) then
440	write (,) 'forcing_GCM2SCM not yet implemented'
441	stop 'in initialization'
442	endif ! forcing_GCM2SCM
443
444	print*,'mxcalc=',mxcalc
445	print*,'zlay=',zlay(mxcalc)
446	print*,'play=',play(mxcalc)
447
448	cAl1 pour SST forced, appellé depuis ocean_forced_noice
449	ts_cur = tsurf ! SST used in read_tsurf1d
450	!=====================================================================
451	! Initialisation de la physique :
452	!=====================================================================
453
454	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
455	!
456	! day_step, iphysiq lus dans gcm.def ci-dessus
457	! timestep: calcule ci-dessous from rday et day_step
458	! ngrid=1
459	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
460	! rday: defini dans suphel.F (86400.)
461	! day_ini: lu dans run.def (dayref)
462	! rlat,rlon lus dans lmdz1d.def
463	! airefi,zcufi,zcvfi initialises au debut de ce programme
464	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
465	day_step = float(nsplit_phys)*day_step/float(iphysiq)
466	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
467	timestep =rday/day_step
468	dtime_frcg = timestep
469	!
470	zcufi=airefi
471	zcvfi=airefi
472
473	call iniphysiq(ngrid,llm,rday,day_ini,timestep,
474	. rlat,rlon,airefi,zcufi,zcvfi,ra,rg,rd,rcpd)
475	print*,'apres iniphysiq'
476
477	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
478	co2_ppm= 330.0
479	solaire=1370.0
480
481	! Ecriture du startphy avant le premier appel a la physique.
482	! On le met juste avant pour avoir acces a tous les champs
483	! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
484
485	if (ok_writedem) then
486
487	!--------------------------------------------------------------------------
488	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
489	! need : qsol fder snow qsurf evap rugos agesno ftsoil
490	!--------------------------------------------------------------------------
491
492	type_ocean = "force"
493	run_off_lic_0(1) = restart_runoff
494	call fonte_neige_init(run_off_lic_0)
495
496	fder=0.
497	snsrf(1,:)=0. ! couverture de neige des sous surface
498	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
499	evap=0.
500	frugs(1,:)=rugos ! couverture de neige des sous surface
501	agesno = xagesno
502	tsoil(:,:,:)=tsurf
503	call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,
504	& evap, frugs, agesno, tsoil)
505
506	!------------------ prepare limit conditions for limit.nc -----------------
507	!-- Ocean force
508
509	print*,'avant phyredem'
510	pctsrf(1,:)=0.
511	if (nat_surf.eq.0.) then
512	pctsrf(1,is_oce)=1.
513	pctsrf(1,is_ter)=0.
514	else
515	pctsrf(1,is_oce)=0.
516	pctsrf(1,is_ter)=1.
517	end if
518
519	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf
520	$ ,pctsrf(1,is_oce),pctsrf(1,is_ter)
521
522	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
523	zpic = zpicinp
524	ftsol=tsurf
525	falb1 = albedo
526	falb2 = albedo
527	rugoro=rugos
528	t_ancien(1,:)=temp(:)
529	q_ancien(1,:)=q(:,1)
530	pbl_tke=1.e-8
531
532	rain_fall=0.
533	snow_fall=0.
534	solsw=0.
535	sollw=0.
536	radsol=0.
537	rnebcon=0.
538	ratqs=0.
539	clwcon=0.
540	zmea=0.
541	zstd=0.
542	zsig=0.
543	zgam=0.
544	zval=0.
545	zthe=0.
546	ema_work1=0.
547	ema_work2=0.
548
549	!------------------------------------------------------------------------
550	! Make file containing restart for the physics (startphy.nc)
551	!
552	! NB: List of the variables to be written by phyredem (via put_field):
553	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
554	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
555	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
556	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
557	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
558	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
559	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,ema_work1,ema_work2
560	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
561	!------------------------------------------------------------------------
562	CAl1 =============== restart option ==========================
563	if (.not.restart) then
564	call phyredem ("startphy.nc")
565	else
566	c (desallocations)
567	print*,'callin surf final'
568	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,
569	& evap, frugs, agesno, tsoil)
570	print*,'after surf final'
571	CALL fonte_neige_final(run_off_lic_0)
572	endif
573
574	ok_writedem=.false.
575	print*,'apres phyredem'
576
577	endif ! ok_writedem
578
579	!------------------------------------------------------------------------
580	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
581	! --------------------------------------------------
582	! NB: List of the variables to be written in limit.nc
583	! (by writelim.F, subroutine of 1DUTILS.h):
584	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
585	! phy_fter,phy_foce,phy_flic,phy_fsic)
586	!------------------------------------------------------------------------
587	do i=1,yd
588	phy_nat(i) = nat_surf
589	phy_alb(i) = albedo
590	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
591	phy_rug(i) = rugos
592	phy_fter(i) = pctsrf(1,is_ter)
593	phy_foce(i) = pctsrf(1,is_oce)
594	phy_fsic(i) = pctsrf(1,is_sic)
595	phy_flic(i) = pctsrf(1,is_lic)
596	enddo
597
598	C fabrication de limit.nc
599	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,
600	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
601
602
603	call phys_state_var_end
604	cAl1
605	if (restart) then
606	print*,'call to restart dyn 1d'
607	Call dyn1deta0("start1dyn.nc",
608	& plev,play,phi,phis,presnivs,
609	& u,v,temp,q,omega2)
610
611	print*,'fnday,annee_ref,day_ref,day_ini',
612	& fnday,annee_ref,day_ref,day_ini
613	c** call ymds2ju(annee_ref,mois,day_ini,heure,day)
614	day = day_ini
615	day_end = day_ini + nday
616	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
617
618	! Print out the actual date of the beginning of the simulation :
619	call ju2ymds(daytime, an, mois, jour, heure)
620	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
621
622	day = int(daytime)
623	time=daytime-day
624
625	print,'*** intialised fields from restart1dyn *****'
626	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
627	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
628	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
629	c raz for safety
630	do l=1,llm
631	dq_dyn(l,1) = 0.
632	enddo
633	endif
634	!Al1 ================ end restart =================================
635	IF (ecrit_slab_oc.eq.1) then
636	open(97,file='div_slab.dat',STATUS='UNKNOWN')
637	elseif (ecrit_slab_oc.eq.0) then
638	open(97,file='div_slab.dat',STATUS='OLD')
639	endif
640	!=====================================================================
641	! START OF THE TEMPORAL LOOP :
642	!=====================================================================
643
644	do while(it.le.nint(fnday*day_step))
645
646	if (prt_level.ge.1) then
647	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',
648	. it,day,time,nint(fnday*day_step)
649	print*,'PAS DE TEMPS ',timestep
650	endif
651	!Al1 demande de restartphy.nc
652	if (it.eq.nint(fnday*day_step)) lastcall=.True.
653
654	!---------------------------------------------------------------------
655	! Interpolation of forcings in time and onto model levels
656	!---------------------------------------------------------------------
657
658	#include "1D_interp_cases.h"
659
660	if (forcing_GCM2SCM) then
661	write (,) 'forcing_GCM2SCM not yet implemented'
662	stop 'in time loop'
663	endif ! forcing_GCM2SCM
664
665	!---------------------------------------------------------------------
666	! Geopotential :
667	!---------------------------------------------------------------------
668
669	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
670	do l = 1, llm-1
671	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1))
672	. (play(l)-play(l+1))/(play(l)+play(l+1))
673	enddo
674
675	!---------------------------------------------------------------------
676	! Call physiq :
677	!---------------------------------------------------------------------
678
679	if (prt_level.ge.1) then
680	print *,' avant physiq : -------- day time ',day,time
681	print *,' temp ',temp
682	print *,' u ',u
683	print *,' q ',q
684	print *,' omega2 ',omega2
685	endif
686	! call writefield_phy('u', u,llm)
687
688	call physiq(ngrid,llm,
689	: firstcall,lastcall,
690	: day,time,timestep,
691	: plev,play,phi,phis,presnivs,clesphy0,
692	: u,v,temp,q,omega2,
693	: du_phys,dv_phys,dt_phys,dq,dpsrf,
694	: dudyn,PVteta)
695	! call writefield_phy('u', u,llm)
696
697	firstcall=.false.
698	if (prt_level.ge.1) then
699	print*,'APRES PHYS'
700	print *,' temp ',temp
701	print *,' q ',q
702	print *,' dq ',dq
703	print*,'dq_dyn',dq_dyn
704	print *,' dt_phys ',dt_phys
705	print *,' dpsrf ',dpsrf
706	print *,' dudyn ',dudyn
707	print *,' PVteta',PVteta
708	endif
709	!---------------------------------------------------------------------
710	! Add physical tendencies :
711	!---------------------------------------------------------------------
712
713	fcoriolis=2.sin(rpixlat/180.)*romega
714
715	if (forcing_radconv) then
716	fcoriolis=0.0
717	dt_cooling=0.0
718	d_th_adv=0.0
719	d_q_adv=0.0
720	endif
721
722	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice) then
723	fcoriolis=0.0
724	endif
725	if(forcing_rico) then
726	dt_cooling=0.
727	endif
728
729	print*, 'fcoriolis ', fcoriolis, rlat(1)
730
731	du_age(1:mxcalc)=
732	: fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
733	dv_age(1:mxcalc)=
734	: -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
735	! call writefield_phy('dv_age' ,dv_age,llm)
736	! call writefield_phy('du_age' ,du_age,llm)
737	! call writefield_phy('du_phys' ,du_phys,llm)
738	! call writefield_phy('u_tend' ,u,llm)
739	! call writefield_phy('u_g' ,ug,llm)
740	u(1:mxcalc)=u(1:mxcalc) + timestep*(
741	: du_phys(1:mxcalc)
742	: +du_age(1:mxcalc) )
743	v(1:mxcalc)=v(1:mxcalc) + timestep*(
744	: dv_phys(1:mxcalc)
745	: +dv_age(1:mxcalc) )
746	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
747	: dq(1:mxcalc,:)
748	: +d_q_adv(1:mxcalc,:) )
749
750	if (prt_level.ge.1) then
751	print *,
752	: 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',
753	: temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
754	print*,du_phys
755	print*, v
756	print*, vg
757	endif
758
759	temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
760	. dt_phys(1:mxcalc)
761	. +d_th_adv(1:mxcalc)
762	. +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
763
764	teta=temp(pzero/play)*rkappa
765
766	!---------------------------------------------------------------------
767	! Add large-scale tendencies (advection, etc) :
768	!---------------------------------------------------------------------
769
770	ccc nrlmd
771	ccc tmpvar=teta
772	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
773	ccc
774	ccc teta(1:mxcalc)=tmpvar(1:mxcalc)
775	ccc tmpvar(:)=q(:,1)
776	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
777	ccc q(1:mxcalc,1)=tmpvar(1:mxcalc)
778	ccc tmpvar(:)=q(:,2)
779	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
780	ccc q(1:mxcalc,2)=tmpvar(1:mxcalc)
781
782	!---------------------------------------------------------------------
783	! Air temperature :
784	!---------------------------------------------------------------------
785	if (lastcall) then
786	print*,'Pas de temps final ',it
787	call ju2ymds(daytime, an, mois, jour, heure)
788	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
789	endif
790
791	c incremente day time
792	c print*,'daytime bef',daytime,1./day_step
793	daytime = daytime+1./day_step
794	!Al1dbg
795	day = int(daytime+0.1/day_step)
796	c time = max(daytime-day,0.0)
797	cAl1&jyg: correction de bug
798	ccc time = real(mod(it,day_step))/day_step
799	time = time_ini/24.+real(mod(it,day_step))/day_step
800	c print*,'daytime nxt time',daytime,time
801	it=it+1
802
803	enddo
804
805	!Al1
806	if (ecrit_slab_oc.ne.-1) close(97)
807
808	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
809	! -------------------------------------
810	call dyn1dredem("restart1dyn.nc",
811	: plev,play,phi,phis,presnivs,
812	: u,v,temp,q,omega2)
813
814	CALL abort_gcm ('lmdz1d ','The End ',0)
815
816	end
817
818	#include "1DUTILS.h"
819	#include "1Dconv.h"
820
821

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

Download in other formats:

Original Format