Context Navigation

lmdz1d.F @ 1795

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

Version testing basee sur la r1794

Testing release based on r1794

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

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

Download in other formats:

Original Format