Context Navigation

lmdz1d.F @ 1607

Last change on this file since 1607 was 1607, checked in by lguez, 13 years ago

Import 1D files. Files added to directory "phy1d" were in directories:

lmdz1d_source_20111207/phy1d_source
lmdz1d_source_20111207/phy1d_source_upd

extracted from:

http://www.lmd.jussieu.fr/~jyg/lmdz1d_source_20111207.tar.gz

File size: 28.3 KB

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

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

Download in other formats:

Original Format