Context Navigation

lmdz1d.F @ 1608

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

In "phy1d/*", replaced obsolete calls to "lnblnk" by calls to "trim".

In "lmdz1d.F", use collector module "ioipsl" rather than specific
module "calendar", which is an internal module of IOIPSL. Removed line
containing only "#" (causes compilation error). Bug fix in call to
"init_phys_lmdz".

File size: 28.3 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	else
276	write (,) 'ERROR : unknown forcing_type ', forcing_type
277	stop 'Forcing_type should be 0,1,2,3 or 40'
278	ENDIF
279	print*,"forcing type=",forcing_type
280
281	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
282	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
283	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
284	! through the common sst_forcing.
285
286	type_ts_forcing = 0
287	if (forcing_toga) type_ts_forcing = 1
288
289	!---------------------------------------------------------------------
290	! Definition of the run
291	!---------------------------------------------------------------------
292
293	call conf_gcm( 99, .TRUE. , clesphy0 )
294	c-----------------------------------------------------------------------
295	c Choix du calendrier
296	c -------------------
297
298	c calend = 'earth_365d'
299	if (calend == 'earth_360d') then
300	call ioconf_calendar('360d')
301	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
302	else if (calend == 'earth_365d') then
303	call ioconf_calendar('noleap')
304	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
305	else if (calend == 'earth_366d') then
306	call ioconf_calendar('all_leap')
307	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
308	else if (calend == 'gregorian') then
309	call ioconf_calendar('gregorian') ! not to be used by normal users
310	write(,)'CALENDRIER CHOISI: Gregorien'
311	else
312	write (,) 'ERROR : unknown calendar ', calend
313	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
314	endif
315	c-----------------------------------------------------------------------
316	c
317	cc Date :
318	c La date est supposee donnee sous la forme [annee, numero du jour dans
319	c l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
320	c On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
321	c Le numero du jour est dans "day". L heure est traitee separement.
322	c La date complete est dans "daytime" (l'unite est le jour).
323	fnday=nday
324	c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
325	IF(forcing_type .EQ. 61) fnday=53100./86400.
326	annee_ref = anneeref
327	mois = 1
328	day_ref = dayref
329	heure = 0.
330	itau_dyn = 0
331	itau_phy = 0
332	call ymds2ju(annee_ref,mois,day_ref,heure,day)
333	day_ini = day
334	day_end = day_ini + nday
335	! Convert the initial date of Toga-Coare to Julian day
336	call ymds2ju
337	$ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
338
339	! Convert the initial date of TWPICE to Julian day
340	call ymds2ju
341	$ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
342	$ ,day_ju_ini_twpi)
343
344	! Convert the initial date of Arm_cu to Julian day
345	call ymds2ju
346	$ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
347	$ ,day_ju_ini_armcu)
348
349	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
350	! Print out the actual date of the beginning of the simulation :
351	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
352	print *,' Time of beginning : ',
353	$ year_print, month_print, day_print, sec_print
354
355	!---------------------------------------------------------------------
356	! Initialization of dimensions, geometry and initial state
357	!---------------------------------------------------------------------
358	call init_phys_lmdz(1,1,llm,1,(/1/))
359	call suphel
360	call initcomgeomphy
361	call infotrac_init
362
363	allocate(q(llm,nqtot))
364	allocate(dq(llm,nqtot))
365	allocate(dq_dyn(llm,nqtot))
366	allocate(d_q_adv(llm,nqtot))
367
368	c
369	c No ozone climatology need be read in this pre-initialization
370	c (phys_state_var_init is called again in physiq)
371	read_climoz = 0
372	c
373	call phys_state_var_init(read_climoz)
374
375	if (ngrid.ne.klon) then
376	print*,'stop in inifis'
377	print*,'Probleme de dimensions :'
378	print*,'ngrid = ',ngrid
379	print*,'klon = ',klon
380	stop
381	endif
382	!!!=====================================================================
383	!!! Feedback forcing values for Gateaux differentiation (al1)
384	!!!=====================================================================
385	!!! Surface Planck forcing bracketing call radiation
386	!! surf_Planck = 0.
387	!! surf_Conv = 0.
388	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
389	!!! a mettre dans le lmdz1d.def ou autre
390	!!
391	!!
392	qsol = qsolinp
393	qsurf = fq_sat(tsurf,psurf/100.)
394	rlat=xlat
395	rlon=xlon
396	day1= day_ini
397	time=daytime-day
398	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
399	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
400
401	ok_flux_surf=.false.
402	fsens=-wtsurfrcpdrho(1)
403	flat=-wqsurfrlvttrho(1)
404
405	! Vertical discretization and pressure levels at half and mid levels:
406
407	pa = 5e4
408	!! preff= 1.01325e5
409	preff = psurf
410	call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
411	sig_s=presnivs/preff
412	plev =ap+bp*psurf
413	play = 0.5*(plev(1:llm)+plev(2:llm+1))
414	ccc zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles
415
416	write(,) '***********************'
417	do l = 1, llm
418	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
419	enddo
420	write(,) '***********************'
421
422	c
423	!=====================================================================
424	! EVENTUALLY, READ FORCING DATA :
425	!=====================================================================
426
427	#include "1D_read_forc_cases.h"
428
429	if (forcing_GCM2SCM) then
430	write (,) 'forcing_GCM2SCM not yet implemented'
431	stop 'in initialization'
432	endif ! forcing_GCM2SCM
433
434	print*,'mxcalc=',mxcalc
435	print*,'zlay=',zlay(mxcalc)
436	print*,'play=',play(mxcalc)
437
438	cAl1 pour SST forced, appellé depuis ocean_forced_noice
439	ts_cur = tsurf ! SST used in read_tsurf1d
440	!=====================================================================
441	! Initialisation de la physique :
442	!=====================================================================
443
444	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
445	!
446	! day_step, iphysiq lus dans gcm.def ci-dessus
447	! timestep: calcule ci-dessous from rday et day_step
448	! ngrid=1
449	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
450	! rday: defini dans suphel.F (86400.)
451	! day_ini: lu dans run.def (dayref)
452	! rlat,rlon lus dans lmdz1d.def
453	! airefi,zcufi,zcvfi initialises au debut de ce programme
454	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
455	day_step = float(nsplit_phys)*day_step/float(iphysiq)
456	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
457	timestep =rday/day_step
458	dtime_frcg = timestep
459	!
460	zcufi=airefi
461	zcvfi=airefi
462
463	call iniphysiq(ngrid,llm,rday,day_ini,timestep,
464	. rlat,rlon,airefi,zcufi,zcvfi,ra,rg,rd,rcpd)
465	print*,'apres iniphysiq'
466
467	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
468	co2_ppm= 330.0
469	solaire=1370.0
470
471	! Ecriture du startphy avant le premier appel a la physique.
472	! On le met juste avant pour avoir acces a tous les champs
473	! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
474
475	if (ok_writedem) then
476
477	!--------------------------------------------------------------------------
478	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
479	! need : qsol fder snow qsurf evap rugos agesno ftsoil
480	!--------------------------------------------------------------------------
481
482	type_ocean = "force"
483	run_off_lic_0(1) = restart_runoff
484	call fonte_neige_init(run_off_lic_0)
485
486	fder=0.
487	snsrf(1,:)=0. ! couverture de neige des sous surface
488	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
489	evap=0.
490	frugs(1,:)=rugos ! couverture de neige des sous surface
491	agesno = xagesno
492	tsoil(:,:,:)=tsurf
493	call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,
494	& evap, frugs, agesno, tsoil)
495
496	!------------------ prepare limit conditions for limit.nc -----------------
497	!-- Ocean force
498
499	print*,'avant phyredem'
500	pctsrf(1,:)=0.
501	if (nat_surf.eq.0.) then
502	pctsrf(1,is_oce)=1.
503	pctsrf(1,is_ter)=0.
504	else
505	pctsrf(1,is_oce)=0.
506	pctsrf(1,is_ter)=1.
507	end if
508
509	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf
510	$ ,pctsrf(1,is_oce),pctsrf(1,is_ter)
511
512	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
513	zpic = zpicinp
514	ftsol=tsurf
515	falb1 = albedo
516	falb2 = albedo
517	rugoro=rugos
518	t_ancien(1,:)=temp(:)
519	q_ancien(1,:)=q(:,1)
520	pbl_tke=1.e-8
521
522	rain_fall=0.
523	snow_fall=0.
524	solsw=0.
525	sollw=0.
526	radsol=0.
527	rnebcon=0.
528	ratqs=0.
529	clwcon=0.
530	zmea=0.
531	zstd=0.
532	zsig=0.
533	zgam=0.
534	zval=0.
535	zthe=0.
536	ema_work1=0.
537	ema_work2=0.
538
539	!------------------------------------------------------------------------
540	! Make file containing restart for the physics (startphy.nc)
541	!
542	! NB: List of the variables to be written by phyredem (via put_field):
543	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
544	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
545	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
546	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
547	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
548	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
549	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,ema_work1,ema_work2
550	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
551	!------------------------------------------------------------------------
552	CAl1 =============== restart option ==========================
553	if (.not.restart) then
554	call phyredem ("startphy.nc")
555	else
556	c (desallocations)
557	print*,'callin surf final'
558	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,
559	& evap, frugs, agesno, tsoil)
560	print*,'after surf final'
561	CALL fonte_neige_final(run_off_lic_0)
562	endif
563
564	ok_writedem=.false.
565	print*,'apres phyredem'
566
567	endif ! ok_writedem
568
569	!------------------------------------------------------------------------
570	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
571	! --------------------------------------------------
572	! NB: List of the variables to be written in limit.nc
573	! (by writelim.F, subroutine of 1DUTILS.h):
574	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
575	! phy_fter,phy_foce,phy_flic,phy_fsic)
576	!------------------------------------------------------------------------
577	do i=1,yd
578	phy_nat(i) = nat_surf
579	phy_alb(i) = albedo
580	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
581	phy_rug(i) = rugos
582	phy_fter(i) = pctsrf(1,is_ter)
583	phy_foce(i) = pctsrf(1,is_oce)
584	phy_fsic(i) = pctsrf(1,is_sic)
585	phy_flic(i) = pctsrf(1,is_lic)
586	enddo
587
588	C fabrication de limit.nc
589	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,
590	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
591
592
593	call phys_state_var_end
594	cAl1
595	if (restart) then
596	print*,'call to restart dyn 1d'
597	Call dyn1deta0("start1dyn.nc",
598	& plev,play,phi,phis,presnivs,
599	& u,v,temp,q,omega2)
600
601	print*,'fnday,annee_ref,day_ref,day_ini',
602	& fnday,annee_ref,day_ref,day_ini
603	c** call ymds2ju(annee_ref,mois,day_ini,heure,day)
604	day = day_ini
605	day_end = day_ini + nday
606	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
607
608	! Print out the actual date of the beginning of the simulation :
609	call ju2ymds(daytime, an, mois, jour, heure)
610	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
611
612	day = int(daytime)
613	time=daytime-day
614
615	print,'*** intialised fields from restart1dyn *****'
616	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
617	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
618	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
619	c raz for safety
620	do l=1,llm
621	dq_dyn(l,1) = 0.
622	enddo
623	endif
624	!Al1 ================ end restart =================================
625	IF (ecrit_slab_oc.eq.1) then
626	open(97,file='div_slab.dat',STATUS='UNKNOWN')
627	elseif (ecrit_slab_oc.eq.0) then
628	open(97,file='div_slab.dat',STATUS='OLD')
629	endif
630	!=====================================================================
631	! START OF THE TEMPORAL LOOP :
632	!=====================================================================
633
634	do while(it.le.nint(fnday*day_step))
635
636	if (prt_level.ge.1) then
637	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',
638	. it,day,time,nint(fnday*day_step)
639	print*,'PAS DE TEMPS ',timestep
640	endif
641	!Al1 demande de restartphy.nc
642	if (it.eq.nint(fnday*day_step)) lastcall=.True.
643
644	!---------------------------------------------------------------------
645	! Interpolation of forcings in time and onto model levels
646	!---------------------------------------------------------------------
647
648	#include "1D_interp_cases.h"
649
650	if (forcing_GCM2SCM) then
651	write (,) 'forcing_GCM2SCM not yet implemented'
652	stop 'in time loop'
653	endif ! forcing_GCM2SCM
654
655	!---------------------------------------------------------------------
656	! Geopotential :
657	!---------------------------------------------------------------------
658
659	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
660	do l = 1, llm-1
661	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1))
662	. (play(l)-play(l+1))/(play(l)+play(l+1))
663	enddo
664
665	!---------------------------------------------------------------------
666	! Call physiq :
667	!---------------------------------------------------------------------
668
669	if (prt_level.ge.1) then
670	print *,' avant physiq : -------- day time ',day,time
671	print *,' temp ',temp
672	print *,' u ',u
673	print *,' q ',q
674	print *,' omega2 ',omega2
675	endif
676	! call writefield_phy('u', u,llm)
677
678	call physiq(ngrid,llm,
679	: firstcall,lastcall,
680	: day,time,timestep,
681	: plev,play,phi,phis,presnivs,clesphy0,
682	: u,v,temp,q,omega2,
683	: du_phys,dv_phys,dt_phys,dq,dpsrf,
684	: dudyn,PVteta)
685	! call writefield_phy('u', u,llm)
686
687	firstcall=.false.
688	if (prt_level.ge.1) then
689	print*,'APRES PHYS'
690	print *,' temp ',temp
691	print *,' q ',q
692	print *,' dq ',dq
693	print*,'dq_dyn',dq_dyn
694	print *,' dt_phys ',dt_phys
695	print *,' dpsrf ',dpsrf
696	print *,' dudyn ',dudyn
697	print *,' PVteta',PVteta
698	endif
699	!---------------------------------------------------------------------
700	! Add physical tendencies :
701	!---------------------------------------------------------------------
702
703	fcoriolis=2.sin(rpirlat(1)/180.)*romega
704
705	if (forcing_radconv) then
706	fcoriolis=0.0
707	dt_cooling=0.0
708	d_th_adv=0.0
709	d_q_adv=0.0
710	endif
711
712	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice) then
713	fcoriolis=0.0
714	endif
715	if(forcing_rico) then
716	dt_cooling=0.
717	endif
718
719	print*, 'fcoriolis ', fcoriolis, rlat(1)
720
721	du_age(1:mxcalc)=
722	: fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
723	dv_age(1:mxcalc)=
724	: -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
725	! call writefield_phy('dv_age' ,dv_age,llm)
726	! call writefield_phy('du_age' ,du_age,llm)
727	! call writefield_phy('du_phys' ,du_phys,llm)
728	! call writefield_phy('u_tend' ,u,llm)
729	! call writefield_phy('u_g' ,ug,llm)
730	u(1:mxcalc)=u(1:mxcalc) + timestep*(
731	: du_phys(1:mxcalc)
732	: +du_age(1:mxcalc) )
733	v(1:mxcalc)=v(1:mxcalc) + timestep*(
734	: dv_phys(1:mxcalc)
735	: +dv_age(1:mxcalc) )
736	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
737	: dq(1:mxcalc,:)
738	: +d_q_adv(1:mxcalc,:) )
739
740	if (prt_level.ge.1) then
741	print *,
742	: 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',
743	: temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
744	print*,du_phys
745	print*, v
746	print*, vg
747	endif
748
749	temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
750	. dt_phys(1:mxcalc)
751	. +d_th_adv(1:mxcalc)
752	. +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
753
754	teta=temp(pzero/play)*rkappa
755
756	!---------------------------------------------------------------------
757	! Add large-scale tendencies (advection, etc) :
758	!---------------------------------------------------------------------
759
760	ccc nrlmd
761	ccc tmpvar=teta
762	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
763	ccc
764	ccc teta(1:mxcalc)=tmpvar(1:mxcalc)
765	ccc tmpvar(:)=q(:,1)
766	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
767	ccc q(1:mxcalc,1)=tmpvar(1:mxcalc)
768	ccc tmpvar(:)=q(:,2)
769	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
770	ccc q(1:mxcalc,2)=tmpvar(1:mxcalc)
771
772	!---------------------------------------------------------------------
773	! Air temperature :
774	!---------------------------------------------------------------------
775	if (lastcall) then
776	print*,'Pas de temps final ',it
777	call ju2ymds(daytime, an, mois, jour, heure)
778	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
779	endif
780
781	c incremente day time
782	c print*,'daytime bef',daytime,1./day_step
783	daytime = daytime+1./day_step
784	!Al1dbg
785	day = int(daytime+0.1/day_step)
786	c time = max(daytime-day,0.0)
787	cAl1&jyg: correction de bug
788	ccc time = real(mod(it,day_step))/day_step
789	time = time_ini/24.+real(mod(it,day_step))/day_step
790	c print*,'daytime nxt time',daytime,time
791	it=it+1
792
793	enddo
794
795	!Al1
796	if (ecrit_slab_oc.ne.-1) close(97)
797
798	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
799	! -------------------------------------
800	call dyn1dredem("restart1dyn.nc",
801	: plev,play,phi,phis,presnivs,
802	: u,v,temp,q,omega2)
803
804	CALL abort_gcm ('lmdz1d ','The End ',0)
805
806	end
807
808	#include "1DUTILS.h"
809	#include "1Dconv.h"
810
811

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

Download in other formats:

Original Format