Context Navigation

lmdz1d.F @ 4236

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

La distinction avec ou sans writelim n'a plus lieu d'être
MP Lefebvre

the writelim/nowritelim distinction is of no further use
MP Lefebvre

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

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

Download in other formats:

Original Format