Context Navigation

lmdz1d.F @ 1780

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

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL

AMMA and Sandu cases included in 1D configuration
MPL

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

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

Download in other formats:

Original Format