Context Navigation

lmdz1d.F @ 1785

Last change on this file since 1785 was 1785, checked in by Ehouarn Millour, 11 years ago

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG

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

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

Download in other formats:

Original Format