Context Navigation

lmdz1d.F @ 1838

Last change on this file since 1838 was 1827, checked in by lguez, 11 years ago

Changed names of variables ema_work1 and ema_work2 to more meaningful
sig1 and w01. Same change in (re)startphy.nc. phyetat0 tries to find
old names ema_work1 and ema_work2 if new names sig1 and w01 are not
found, so the program can run with an old restartphy.nc. restartphy.nc
is modified compared to the previous SVN revision because of the change of
names but the data content is not modified (this can be checked with
max_diff_nc.sh -i).

File size: 34.7 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	sig1=0.
635	w01=0.
636	u_ancien(1,:)=u(:)
637	v_ancien(1,:)=v(:)
638
639	!------------------------------------------------------------------------
640	! Make file containing restart for the physics (startphy.nc)
641	!
642	! NB: List of the variables to be written by phyredem (via put_field):
643	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
644	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
645	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
646	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
647	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
648	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
649	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,sig1,w01
650	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
651	!------------------------------------------------------------------------
652	CAl1 =============== restart option ==========================
653	if (.not.restart) then
654	call phyredem ("startphy.nc")
655	else
656	c (desallocations)
657	print*,'callin surf final'
658	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,
659	& evap, frugs, agesno, tsoil)
660	print*,'after surf final'
661	CALL fonte_neige_final(run_off_lic_0)
662	endif
663
664	ok_writedem=.false.
665	print*,'apres phyredem'
666
667	endif ! ok_writedem
668
669	!------------------------------------------------------------------------
670	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
671	! --------------------------------------------------
672	! NB: List of the variables to be written in limit.nc
673	! (by writelim.F, subroutine of 1DUTILS.h):
674	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
675	! phy_fter,phy_foce,phy_flic,phy_fsic)
676	!------------------------------------------------------------------------
677	do i=1,yd
678	phy_nat(i) = nat_surf
679	phy_alb(i) = albedo
680	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
681	phy_rug(i) = rugos
682	phy_fter(i) = pctsrf(1,is_ter)
683	phy_foce(i) = pctsrf(1,is_oce)
684	phy_fsic(i) = pctsrf(1,is_sic)
685	phy_flic(i) = pctsrf(1,is_lic)
686	enddo
687
688	C fabrication de limit.nc
689	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,
690	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
691
692
693	call phys_state_var_end
694	cAl1
695	if (restart) then
696	print*,'call to restart dyn 1d'
697	Call dyn1deta0("start1dyn.nc",
698	& plev,play,phi,phis,presnivs,
699	& u,v,temp,q,omega2)
700
701	print*,'fnday,annee_ref,day_ref,day_ini',
702	& fnday,annee_ref,day_ref,day_ini
703	c** call ymds2ju(annee_ref,mois,day_ini,heure,day)
704	day = day_ini
705	day_end = day_ini + nday
706	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
707
708	! Print out the actual date of the beginning of the simulation :
709	call ju2ymds(daytime, an, mois, jour, heure)
710	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
711
712	day = int(daytime)
713	time=daytime-day
714
715	print,'*** intialised fields from restart1dyn *****'
716	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
717	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
718	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
719	c raz for safety
720	do l=1,llm
721	dq_dyn(l,1) = 0.
722	enddo
723	endif
724	!Al1 ================ end restart =================================
725	IF (ecrit_slab_oc.eq.1) then
726	open(97,file='div_slab.dat',STATUS='UNKNOWN')
727	elseif (ecrit_slab_oc.eq.0) then
728	open(97,file='div_slab.dat',STATUS='OLD')
729	endif
730	!=====================================================================
731	! START OF THE TEMPORAL LOOP :
732	!=====================================================================
733
734	do while(it.le.nint(fnday*day_step))
735
736	if (prt_level.ge.1) then
737	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',
738	. it,day,time,nint(fnday*day_step)
739	print*,'PAS DE TEMPS ',timestep
740	endif
741	!Al1 demande de restartphy.nc
742	if (it.eq.nint(fnday*day_step)) lastcall=.True.
743
744	!---------------------------------------------------------------------
745	! Interpolation of forcings in time and onto model levels
746	!---------------------------------------------------------------------
747
748	#include "1D_interp_cases.h"
749
750	if (forcing_GCM2SCM) then
751	write (,) 'forcing_GCM2SCM not yet implemented'
752	stop 'in time loop'
753	endif ! forcing_GCM2SCM
754
755	!---------------------------------------------------------------------
756	! Geopotential :
757	!---------------------------------------------------------------------
758
759	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
760	do l = 1, llm-1
761	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1))
762	. (play(l)-play(l+1))/(play(l)+play(l+1))
763	enddo
764
765	!---------------------------------------------------------------------
766	! Listing output for debug prt_level>=1
767	!---------------------------------------------------------------------
768	if (prt_level>=1) then
769	print *,' avant physiq : -------- day time ',day,time
770	write(,) 'firstcall,lastcall,phis',
771	: firstcall,lastcall,phis
772	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',
773	: 'presniv','plev','play','phi'
774	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,
775	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
776	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',
777	: 'presniv','u','v','temp','q1','q2','omega2'
778	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,
779	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
780	endif
781
782	!---------------------------------------------------------------------
783	! Call physiq :
784	!---------------------------------------------------------------------
785
786	call physiq(ngrid,llm,
787	: firstcall,lastcall,
788	: day,time,timestep,
789	: plev,play,phi,phis,presnivs,clesphy0,
790	: u,v,temp,q,omega2,
791	: du_phys,dv_phys,dt_phys,dq,dpsrf,
792	: dudyn,PVteta)
793	firstcall=.false.
794
795	!---------------------------------------------------------------------
796	! Listing output for debug prt_level>=1
797	!---------------------------------------------------------------------
798	if (prt_level>=1) then
799	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',
800	: 'presniv','plev','play','phi'
801	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,
802	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
803	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',
804	: 'presniv','u','v','temp','q1','q2','omega2'
805	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,
806	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
807	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',
808	: 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
809	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,
810	: presnivs(l),86400du_phys(l),86400dv_phys(l),
811	: 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
812	write(,) 'dpsrf',dpsrf
813	endif
814	!---------------------------------------------------------------------
815	! Add physical tendencies :
816	!---------------------------------------------------------------------
817
818	fcoriolis=2.sin(rpixlat/180.)*romega
819	if (forcing_radconv) then
820	fcoriolis=0.0
821	dt_cooling=0.0
822	d_th_adv=0.0
823	d_q_adv=0.0
824	endif
825
826	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice
827	: .or.forcing_amma) then
828	fcoriolis=0.0 ; ug=0. ; vg=0.
829	endif
830	if(forcing_rico) then
831	dt_cooling=0.
832	endif
833
834	print*, 'fcoriolis ', fcoriolis, xlat
835
836	du_age(1:mxcalc)=
837	: fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
838	dv_age(1:mxcalc)=
839	: -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
840
841	!!!!!!!!!!!!!!!!!!!!!!!!
842	! Geostrophic wind
843	!!!!!!!!!!!!!!!!!!!!!!!!
844	sfdt = sin(0.5fcoriolistimestep)
845	cfdt = cos(0.5fcoriolistimestep)
846	!
847	du_age(1:mxcalc)= -2.sfdt/timestep
848	: (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -
849	: cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
850	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
851	!
852	dv_age(1:mxcalc)= -2.sfdt/timestep
853	: (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +
854	: sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
855	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
856	!
857	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
858	! call writefield_phy('dv_age' ,dv_age,llm)
859	! call writefield_phy('du_age' ,du_age,llm)
860	! call writefield_phy('du_phys' ,du_phys,llm)
861	! call writefield_phy('u_tend' ,u,llm)
862	! call writefield_phy('u_g' ,ug,llm)
863	!
864	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
865	!! Increment state variables
866	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
867	u(1:mxcalc)=u(1:mxcalc) + timestep*(
868	: du_phys(1:mxcalc)
869	: +du_age(1:mxcalc) )
870	v(1:mxcalc)=v(1:mxcalc) + timestep*(
871	: dv_phys(1:mxcalc)
872	: +dv_age(1:mxcalc) )
873	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
874	: dq(1:mxcalc,:)
875	: +d_q_adv(1:mxcalc,:) )
876
877	if (prt_level.ge.1) then
878	print *,
879	: 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',
880	: temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
881	print*,du_phys
882	print*, v
883	print*, vg
884	endif
885
886	temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
887	. dt_phys(1:mxcalc)
888	. +d_th_adv(1:mxcalc)
889	. +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
890
891	teta=temp(pzero/play)*rkappa
892	!
893	!---------------------------------------------------------------------
894	! Nudge soil temperature if requested
895	!---------------------------------------------------------------------
896
897	IF (nudge_tsoil) THEN
898	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)
899	. -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
900	ENDIF
901
902	!---------------------------------------------------------------------
903	! Add large-scale tendencies (advection, etc) :
904	!---------------------------------------------------------------------
905
906	ccc nrlmd
907	ccc tmpvar=teta
908	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
909	ccc
910	ccc teta(1:mxcalc)=tmpvar(1:mxcalc)
911	ccc tmpvar(:)=q(:,1)
912	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
913	ccc q(1:mxcalc,1)=tmpvar(1:mxcalc)
914	ccc tmpvar(:)=q(:,2)
915	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
916	ccc q(1:mxcalc,2)=tmpvar(1:mxcalc)
917
918	!---------------------------------------------------------------------
919	! Air temperature :
920	!---------------------------------------------------------------------
921	if (lastcall) then
922	print*,'Pas de temps final ',it
923	call ju2ymds(daytime, an, mois, jour, heure)
924	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
925	endif
926
927	c incremente day time
928	c print*,'daytime bef',daytime,1./day_step
929	daytime = daytime+1./day_step
930	!Al1dbg
931	day = int(daytime+0.1/day_step)
932	c time = max(daytime-day,0.0)
933	cAl1&jyg: correction de bug
934	ccc time = real(mod(it,day_step))/day_step
935	time = time_ini/24.+real(mod(it,day_step))/day_step
936	c print*,'daytime nxt time',daytime,time
937	it=it+1
938
939	enddo
940
941	!Al1
942	if (ecrit_slab_oc.ne.-1) close(97)
943
944	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
945	! -------------------------------------
946	call dyn1dredem("restart1dyn.nc",
947	: plev,play,phi,phis,presnivs,
948	: u,v,temp,q,omega2)
949
950	CALL abort_gcm ('lmdz1d ','The End ',0)
951
952	end
953
954	#include "1DUTILS.h"
955	#include "1Dconv.h"
956
957

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

Download in other formats:

Original Format