Context Navigation

scm.f90 @ 5451

Last change on this file since 5451 was 5433, checked in by evignon, 6 days ago

modification et correction de la prise en compte de l'advection verticale dans le 1D.
Il y a 3 aspects:

on tient compte de l'advection verticale des 3 phases de l'eau
on corrige l'application de la tendance qui etait calculee avant l'appel de la physique mais appliquee apres
on utilise desormais un schema amont par defaut (le schema centre, mis sous flag, conduisait a des oscillations). C'est un peu fort. Implementation a venir d'un schema van leer

ATTENTION: on perd la convergence sur les cas 1D avec vitesse verticale grande echelle (RICO, SANDU ...)

File size: 39.9 KB

Line
1	SUBROUTINE scm
2
3	USE flux_arp_mod_h
4	USE compbl_mod_h
5	USE clesphys_mod_h
6	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
7	USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
8	clwcon, detr_therm, &
9	qsol, fevap, z0m, z0h, agesno, &
10	du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
11	falb_dir, falb_dif, &
12	ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
13	rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
14	solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, &
15	wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
16	wake_deltaq, wake_deltat, wake_s, awake_s, wake_dens, &
17	awake_dens, cv_gen, wake_cstar, &
18	zgam, zmax0, zmea, zpic, zsig, &
19	zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
20	ql_ancien, qs_ancien, qbs_ancien, cf_ancien, rvc_ancien, &
21	prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
22	u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
23
24
25	USE dimphy
26	USE surface_data, only : type_ocean,ok_veget
27	USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
28	pbl_surface_final
29	USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
30
31	USE infotrac ! new
32	USE control_mod
33	USE indice_sol_mod
34	USE phyaqua_mod
35	USE mod_1D_cases_read_std
36	USE print_control_mod, ONLY: lunout, prt_level
37	USE iniphysiq_mod, ONLY: iniphysiq
38	USE mod_const_mpi, ONLY: comm_lmdz
39	USE physiq_mod, ONLY: physiq
40	USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
41	preff, aps, bps, pseudoalt, scaleheight
42	USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
43	itau_dyn, itau_phy, start_time, year_len
44	USE phys_cal_mod, ONLY : year_len_phys_cal_mod => year_len
45	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_OUTPUTPHYSSCM
46
47	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
48	USE dimsoil_mod_h, ONLY: nsoilmx
49	USE yomcst_mod_h
50	USE tsoilnudge_mod_h
51	USE fcg_gcssold_mod_h
52	USE compar1d_mod_h
53	USE date_cas_mod_h
54	implicit none
55
56	!=====================================================================
57	! DECLARATIONS
58	!=====================================================================
59
60	!---------------------------------------------------------------------
61	! Externals
62	!---------------------------------------------------------------------
63	external fq_sat
64	real fq_sat
65
66	!---------------------------------------------------------------------
67	! Arguments d' initialisations de la physique (USER DEFINE)
68	!---------------------------------------------------------------------
69
70	integer, parameter :: ngrid=1
71	real :: zcufi = 1.
72	real :: zcvfi = 1.
73	real :: fnday
74	real :: day, daytime
75	real :: day1
76	real :: heure
77	integer :: jour
78	integer :: mois
79	integer :: an
80
81	!---------------------------------------------------------------------
82	! Declarations related to forcing and initial profiles
83	!---------------------------------------------------------------------
84
85	integer :: kmax = llm
86	integer llm700,nq1,nq2,iflag_1d_vert_adv
87	INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
88	real timestep, frac
89	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
90	real uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
91	real ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
92	real dqtdxls(nlev_max),dqtdyls(nlev_max)
93	real dqtdtls(nlev_max),thlpcar(nlev_max)
94	real qprof(nlev_max,nqmx)
95
96	! integer :: forcing_type
97	logical :: forcing_les = .false.
98	logical :: forcing_armcu = .false.
99	logical :: forcing_rico = .false.
100	logical :: forcing_radconv = .false.
101	logical :: forcing_toga = .false.
102	logical :: forcing_twpice = .false.
103	logical :: forcing_amma = .false.
104	logical :: forcing_dice = .false.
105	logical :: forcing_gabls4 = .false.
106
107	logical :: forcing_GCM2SCM = .false.
108	logical :: forcing_GCSSold = .false.
109	logical :: forcing_sandu = .false.
110	logical :: forcing_astex = .false.
111	logical :: forcing_fire = .false.
112	logical :: forcing_case = .false.
113	logical :: forcing_case2 = .false.
114	logical :: forcing_SCM = .false.
115
116	!flag forcings
117	logical :: nudge_wind=.true.
118	logical :: nudge_thermo=.false.
119	logical :: cptadvw=.true.
120
121
122	!=====================================================================
123	! DECLARATIONS FOR EACH CASE
124	!=====================================================================
125	!
126	INCLUDE "1D_decl_cases.h"
127	!
128	!---------------------------------------------------------------------
129	! Declarations related to nudging
130	!---------------------------------------------------------------------
131	integer :: nudge_max
132	parameter (nudge_max=9)
133	integer :: inudge_RHT=1
134	integer :: inudge_UV=2
135	logical :: nudge(nudge_max)
136	real :: t_targ(llm)
137	real :: rh_targ(llm)
138	real :: u_targ(llm)
139	real :: v_targ(llm)
140	!
141	!---------------------------------------------------------------------
142	! Declarations related to vertical discretization:
143	!---------------------------------------------------------------------
144	real :: pzero=1.e5
145	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
146	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
147
148	!---------------------------------------------------------------------
149	! Declarations related to variables
150	!---------------------------------------------------------------------
151
152	real :: phi(llm)
153	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
154	REAL rot(1, llm) ! relative vorticity, in s-1
155	real :: rlat_rad(1),rlon_rad(1)
156	real :: omega(llm),omega2(llm),rho(llm+1)
157	real :: ug(llm),vg(llm),fcoriolis
158	real :: sfdt, cfdt
159	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
160	real :: w_adv(llm),z_adv(llm)
161	real :: d_t_vert_adv(llm),d_u_vert_adv(llm),d_v_vert_adv(llm)
162	real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
163	real :: d_u_nudge(llm),d_v_nudge(llm)
164	real :: d_u_age(llm),d_v_age(llm)
165	real :: alpha
166	real :: ttt
167
168	REAL, ALLOCATABLE, DIMENSION(:,:):: q
169	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
170	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_vert_adv
171	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
172	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
173
174	!---------------------------------------------------------------------
175	! Initialization of surface variables
176	!---------------------------------------------------------------------
177	real :: run_off_lic_0(1)
178	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
179	real :: tsoil(1,nsoilmx,nbsrf)
180
181	!---------------------------------------------------------------------
182	! Call to phyredem
183	!---------------------------------------------------------------------
184	logical :: ok_writedem =.true.
185	real :: sollw_in = 0.
186	real :: solsw_in = 0.
187
188	!---------------------------------------------------------------------
189	! Call to physiq
190	!---------------------------------------------------------------------
191	logical :: firstcall=.true.
192	logical :: lastcall=.false.
193	real :: phis(1) = 0.0
194	real :: dpsrf(1)
195
196	!---------------------------------------------------------------------
197	! Initializations of boundary conditions
198	!---------------------------------------------------------------------
199	real, allocatable :: phy_nat (:) ! 0=ocean libre,1=land,2=glacier,3=banquise
200	real, allocatable :: phy_alb (:) ! Albedo land only (old value condsurf_jyg=0.3)
201	real, allocatable :: phy_sst (:) ! SST (will not be used; cf read_tsurf1d.F)
202	real, allocatable :: phy_bil (:) ! Ne sert que pour les slab_ocean
203	real, allocatable :: phy_rug (:) ! Longueur rugosite utilisee sur land only
204	real, allocatable :: phy_ice (:) ! Fraction de glace
205	real, allocatable :: phy_fter(:) ! Fraction de terre
206	real, allocatable :: phy_foce(:) ! Fraction de ocean
207	real, allocatable :: phy_fsic(:) ! Fraction de glace
208	real, allocatable :: phy_flic(:) ! Fraction de glace
209
210	!---------------------------------------------------------------------
211	! Fichiers et d'autres variables
212	!---------------------------------------------------------------------
213	integer :: k,l,i,it=1,mxcalc
214	integer :: nsrf
215	integer jcode
216	INTEGER read_climoz
217	integer :: it_end ! iteration number of the last call
218	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
219	data ecrit_slab_oc/-1/
220
221	! if flag_inhib_forcing = 0, tendencies of forcing are added
222	! <> 0, tendencies of forcing are not added
223	INTEGER :: flag_inhib_forcing = 0
224
225
226	print*,'VOUS ENTREZ DANS LE 1D FORMAT STANDARD'
227
228	!=====================================================================
229	! INITIALIZATIONS
230	!=====================================================================
231	du_phys(:)=0.
232	dv_phys(:)=0.
233	dt_phys(:)=0.
234	d_t_vert_adv(:)=0.
235	d_u_vert_adv(:)=0.
236	d_v_vert_adv(:)=0.
237	dt_cooling(:)=0.
238	d_t_adv(:)=0.
239	d_t_nudge(:)=0.
240	d_u_nudge(:)=0.
241	d_v_nudge(:)=0.
242	d_u_adv(:)=0.
243	d_v_adv(:)=0.
244	d_u_age(:)=0.
245	d_v_age(:)=0.
246
247
248	! Initialization of Common turb_forcing
249	dtime_frcg = 0.
250	Turb_fcg_gcssold=.false.
251	hthturb_gcssold = 0.
252	hqturb_gcssold = 0.
253
254
255
256
257	!---------------------------------------------------------------------
258	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
259	!---------------------------------------------------------------------
260	call conf_unicol
261	!Al1 moves this gcssold var from common fcg_gcssold to
262	Turb_fcg_gcssold = xTurb_fcg_gcssold
263	! --------------------------------------------------------------------
264	close(1)
265	write(,) 'lmdz1d.def lu => unicol.def'
266
267	forcing_SCM = .true.
268	year_ini_cas=1997
269	! It is possible that those parameters are run twice.
270	! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
271
272
273	call getin('anneeref',year_ini_cas)
274	call getin('dayref',day_deb)
275	mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
276	call getin('time_ini',heure_ini_cas)
277	print*,'NATURE DE LA SURFACE ',nat_surf
278
279	! Initialization of the logical switch for nudging
280
281	jcode = iflag_nudge
282	do i = 1,nudge_max
283	nudge(i) = mod(jcode,10) .ge. 1
284	jcode = jcode/10
285	enddo
286	! numerical scheme for vertical advection
287
288	iflag_1d_vert_adv=1
289	call getin('iflag_1d_vert_adv',iflag_1d_vert_adv)
290
291	!-----------------------------------------------------------------------
292	! Definition of the run
293	!-----------------------------------------------------------------------
294
295	call conf_gcm( 99, .TRUE. )
296
297	!-----------------------------------------------------------------------
298	allocate( phy_nat (year_len)) ! 0=ocean libre,1=land,2=glacier,3=banquise
299	phy_nat(:)=0.0
300	allocate( phy_alb (year_len)) ! Albedo land only (old value condsurf_jyg=0.3)
301	allocate( phy_sst (year_len)) ! SST (will not be used; cf read_tsurf1d.F)
302	allocate( phy_bil (year_len)) ! Ne sert que pour les slab_ocean
303	phy_bil(:)=1.0
304	allocate( phy_rug (year_len)) ! Longueur rugosite utilisee sur land only
305	allocate( phy_ice (year_len)) ! Fraction de glace
306	phy_ice(:)=0.0
307	allocate( phy_fter(year_len)) ! Fraction de terre
308	phy_fter(:)=0.0
309	allocate( phy_foce(year_len)) ! Fraction de ocean
310	phy_foce(:)=0.0
311	allocate( phy_fsic(year_len)) ! Fraction de glace
312	phy_fsic(:)=0.0
313	allocate( phy_flic(year_len)) ! Fraction de glace
314	phy_flic(:)=0.0
315
316
317	!-----------------------------------------------------------------------
318	! Choix du calendrier
319	! -------------------
320
321	! calend = 'earth_365d'
322	if (calend == 'earth_360d') then
323	call ioconf_calendar('360_day')
324	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
325	else if (calend == 'earth_365d') then
326	call ioconf_calendar('noleap')
327	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
328	else if (calend == 'earth_366d') then
329	call ioconf_calendar('all_leap')
330	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
331	else if (calend == 'gregorian') then
332	stop 'gregorian calend should not be used by normal user'
333	call ioconf_calendar('gregorian') ! not to be used by normal users
334	write(,)'CALENDRIER CHOISI: Gregorien'
335	else
336	write (,) 'ERROR : unknown calendar ', calend
337	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
338	endif
339	!-----------------------------------------------------------------------
340	!
341	!c Date :
342	! La date est supposee donnee sous la forme [annee, numero du jour dans
343	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
344	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
345	! Le numero du jour est dans "day". L heure est traitee separement.
346	! La date complete est dans "daytime" (l'unite est le jour).
347
348
349	if (nday>0) then
350	fnday=nday
351	else
352	fnday=-nday/float(day_step)
353	endif
354	print *,'fnday=',fnday
355	! start_time doit etre en FRACTION DE JOUR
356	start_time=time_ini/24.
357
358	annee_ref = anneeref
359	mois = 1
360	day_ref = dayref
361	heure = 0.
362	itau_dyn = 0
363	itau_phy = 0
364	call ymds2ju(annee_ref,mois,day_ref,heure,day)
365	day_ini = int(day)
366	day_end = day_ini + int(fnday)
367
368	! Convert the initial date to Julian day
369	day_ini_cas=day_deb
370	print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas
371	call ymds2ju &
372	& (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600 &
373	& ,day_ju_ini_cas)
374	print*,'time case 2',day_ini_cas,day_ju_ini_cas
375	daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation
376
377	! Print out the actual date of the beginning of the simulation :
378	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
379	print *,' Time of beginning : ', &
380	& year_print, month_print, day_print, sec_print
381
382	!---------------------------------------------------------------------
383	! Initialization of dimensions, geometry and initial state
384	!---------------------------------------------------------------------
385	! call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
386	! but we still need to initialize dimphy module (klon,klev,etc.) here.
387	call init_dimphy1D(1,llm)
388	call suphel
389	call init_infotrac
390
391	if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
392	allocate(q(llm,nqtot)) ; q(:,:)=0.
393	allocate(dq(llm,nqtot))
394	allocate(d_q_vert_adv(llm,nqtot))
395	allocate(d_q_adv(llm,nqtot))
396	allocate(d_q_nudge(llm,nqtot))
397
398	q(:,:) = 0.
399	dq(:,:) = 0.
400	d_q_vert_adv(:,:) = 0.
401	d_q_adv(:,:) = 0.
402	d_q_nudge(:,:) = 0.
403
404	! No ozone climatology need be read in this pre-initialization
405	! (phys_state_var_init is called again in physiq)
406	read_climoz = 0
407	nsw=6
408
409	call phys_state_var_init(read_climoz)
410
411	if (ngrid.ne.klon) then
412	print*,'stop in inifis'
413	print*,'Probleme de dimensions :'
414	print*,'ngrid = ',ngrid
415	print*,'klon = ',klon
416	stop
417	endif
418	!!!=====================================================================
419	!!! Feedback forcing values for Gateaux differentiation (al1)
420	!!!=====================================================================
421
422	qsol = qsolinp
423	qsurf = fq_sat(tsurf,psurf/100.)
424	beta_aridity(:,:) = beta_surf
425	day1= day_ini
426	time=daytime-day
427	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
428	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
429
430	!! mpl et jyg le 22/08/2012 :
431	!! pour que les cas a flux de surface imposes marchent
432	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
433	fsens=-wtsurfrcpdrho(1)
434	flat=-wqsurfrlvttrho(1)
435	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
436	ENDIF
437	print*,'Flux sol ',fsens,flat
438
439	! Vertical discretization and pressure levels at half and mid levels:
440
441	pa = 5e4
442	!! preff= 1.01325e5
443	preff = psurf
444	IF (ok_old_disvert) THEN
445	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
446	print *,'On utilise disvert0'
447	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
448	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
449	scaleheight=8.
450	pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
451	ELSE
452	call disvert()
453	print *,'On utilise disvert'
454	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
455	! Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
456	ENDIF
457
458	sig_s=presnivs/preff
459	plev =ap+bp*psurf
460	play = 0.5*(plev(1:llm)+plev(2:llm+1))
461	zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles.
462
463	IF (forcing_type .eq. 59) THEN
464	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
465	write(,) '***********************'
466	do l = 1, llm
467	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
468	if (trouve_700 .and. play(l).le.70000) then
469	llm700=l
470	print *,'llm700,play=',llm700,play(l)/100.
471	trouve_700= .false.
472	endif
473	enddo
474	write(,) '***********************'
475	ENDIF
476
477	!
478	!=====================================================================
479	! EVENTUALLY, READ FORCING DATA :
480	!=====================================================================
481
482	INCLUDE "1D_read_forc_cases.h"
483
484	if (forcing_GCM2SCM) then
485	write (,) 'forcing_GCM2SCM not yet implemented'
486	stop 'in initialization'
487	endif ! forcing_GCM2SCM
488
489
490	!=====================================================================
491	! Initialisation de la physique :
492	!=====================================================================
493
494	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
495	!
496	! day_step, iphysiq lus dans gcm.def ci-dessus
497	! timestep: calcule ci-dessous from rday et day_step
498	! ngrid=1
499	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
500	! rday: defini dans suphel.F (86400.)
501	! day_ini: lu dans run.def (dayref)
502	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
503	! airefi,zcufi,zcvfi initialises au debut de ce programme
504	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
505
506
507	day_step = float(nsplit_phys)*day_step/float(iphysiq)
508	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
509	timestep =rday/day_step
510	dtime_frcg = timestep
511	!
512	zcufi=airefi
513	zcvfi=airefi
514	!
515	rlat_rad(1)=xlat*rpi/180.
516	rlon_rad(1)=xlon*rpi/180.
517
518	! iniphysiq will call iniaqua who needs year_len from phys_cal_mod
519	year_len_phys_cal_mod=year_len
520
521	! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
522	! e.g. for cell boundaries, which are meaningless in 1D; so pad these
523	! with '0.' when necessary
524
525	call iniphysiq(iim,jjm,llm, &
526	1,comm_lmdz, &
527	rday,day_ini,timestep, &
528	(/rlat_rad(1),0./),(/0./), &
529	(/0.,0./),(/rlon_rad(1),0./), &
530	(/ (/airefi,0./),(/0.,0./) /), &
531	(/zcufi,0.,0.,0./), &
532	(/zcvfi,0./), &
533	ra,rg,rd,rcpd,1)
534	print*,'apres iniphysiq'
535
536	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
537	co2_ppm= 330.0
538	solaire=1370.0
539
540	! Ecriture du startphy avant le premier appel a la physique.
541	! On le met juste avant pour avoir acces a tous les champs
542
543	if (ok_writedem) then
544
545	!--------------------------------------------------------------------------
546	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
547	! need : qsol fder snow qsurf evap rugos agesno ftsoil
548	!--------------------------------------------------------------------------
549
550	type_ocean = "force"
551	run_off_lic_0(1) = restart_runoff
552	call fonte_neige_init(run_off_lic_0)
553
554	fder=0.
555	snsrf(1,:)=snowmass ! masse de neige des sous surface
556	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
557	fevap=0.
558	z0m(1,:)=rugos ! couverture de neige des sous surface
559	z0h(1,:)=rugosh ! couverture de neige des sous surface
560	agesno = xagesno
561	tsoil(:,:,:)=tsurf
562	!-----------------------------------------------------------------------
563	call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
564
565	!------------------ prepare limit conditions for limit.nc -----------------
566	!-- Ocean force
567
568	print*,'avant phyredem'
569	pctsrf(1,:)=0.
570	if (nat_surf.eq.0.) then
571	pctsrf(1,is_oce)=1.
572	pctsrf(1,is_ter)=0.
573	pctsrf(1,is_lic)=0.
574	pctsrf(1,is_sic)=0.
575	else if (nat_surf .eq. 1) then
576	pctsrf(1,is_oce)=0.
577	pctsrf(1,is_ter)=1.
578	pctsrf(1,is_lic)=0.
579	pctsrf(1,is_sic)=0.
580	else if (nat_surf .eq. 2) then
581	pctsrf(1,is_oce)=0.
582	pctsrf(1,is_ter)=0.
583	pctsrf(1,is_lic)=1.
584	pctsrf(1,is_sic)=0.
585	else if (nat_surf .eq. 3) then
586	pctsrf(1,is_oce)=0.
587	pctsrf(1,is_ter)=0.
588	pctsrf(1,is_lic)=0.
589	pctsrf(1,is_sic)=1.
590
591	end if
592
593
594	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf &
595	& ,pctsrf(1,is_oce),pctsrf(1,is_ter)
596
597	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
598	zpic = zpicinp
599	ftsol=tsurf
600	falb_dir=albedo
601	falb_dif=albedo
602	rugoro=rugos
603	t_ancien(1,:)=temp(:)
604	q_ancien(1,:)=q(:,1)
605	ql_ancien = 0.
606	qs_ancien = 0.
607	prlw_ancien = 0.
608	prsw_ancien = 0.
609	prw_ancien = 0.
610	IF ( ok_bs ) THEN
611	qbs_ancien = 0.
612	prbsw_ancien = 0.
613	ENDIF
614	IF ( ok_ice_supersat ) THEN
615	cf_ancien = 0.
616	rvc_ancien = 0.
617	ENDIF
618	rain_fall=0.
619	snow_fall=0.
620	solsw=0.
621	solswfdiff=0.
622	sollw=0.
623	sollwdown=rsigmatsurf*4
624	radsol=0.
625	rnebcon=0.
626	ratqs=0.
627	clwcon=0.
628	zmax0 = 0.
629	zmea=zsurf
630	zstd=0.
631	zsig=0.
632	zgam=0.
633	zval=0.
634	zthe=0.
635	sig1=0.
636	w01=0.
637	!
638	wake_deltaq = 0.
639	wake_deltat = 0.
640	wake_delta_pbl_TKE(:,:,:) = 0.
641	delta_tsurf = 0.
642	wake_fip = 0.
643	wake_pe = 0.
644	wake_s = 0.
645	awake_s = 0.
646	wake_dens = 0.
647	awake_dens = 0.
648	cv_gen = 0.
649	wake_cstar = 0.
650	ale_bl = 0.
651	ale_bl_trig = 0.
652	alp_bl = 0.
653	IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
654	IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
655	entr_therm = 0.
656	detr_therm = 0.
657	f0 = 0.
658	fm_therm = 0.
659	u_ancien(1,:)=u(:)
660	v_ancien(1,:)=v(:)
661
662	u10m=0.
663	v10m=0.
664	ale_wake=0.
665	ale_bl_stat=0.
666	ratqs_inter_(:,:)= 0.002
667
668	!------------------------------------------------------------------------
669	! Make file containing restart for the physics (startphy.nc)
670	!
671	! NB: List of the variables to be written by phyredem (via put_field):
672	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
673	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
674	! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
675	! radsol,solsw,solswfdiff,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
676	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
677	! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
678	! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
679	! wake_deltat,wake_deltaq,wake_s,awake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
680	! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
681	!
682	! NB2: The content of the startphy.nc file depends on some flags defined in
683	! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
684	! to be set at some arbitratry convenient values.
685	!------------------------------------------------------------------------
686	!Al1 =============== restart option ======================================
687	iflag_physiq=0
688	call getin('iflag_physiq',iflag_physiq)
689
690	if (.not.restart) then
691	iflag_pbl = 5
692	call phyredem ("startphy.nc")
693	else
694	! (desallocations)
695	print*,'callin surf final'
696	call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
697	print*,'after surf final'
698	CALL fonte_neige_final(run_off_lic_0)
699	endif
700
701	ok_writedem=.false.
702	print*,'apres phyredem'
703
704	endif ! ok_writedem
705
706	!------------------------------------------------------------------------
707	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
708	! --------------------------------------------------
709	! NB: List of the variables to be written in limit.nc
710	! (by writelim.F, subroutine of 1DUTILS.h):
711	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
712	! phy_fter,phy_foce,phy_flic,phy_fsic)
713	!------------------------------------------------------------------------
714	do i=1,year_len
715	phy_nat(i) = nat_surf
716	phy_alb(i) = albedo
717	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
718	phy_rug(i) = rugos
719	phy_fter(i) = pctsrf(1,is_ter)
720	phy_foce(i) = pctsrf(1,is_oce)
721	phy_fsic(i) = pctsrf(1,is_sic)
722	phy_flic(i) = pctsrf(1,is_lic)
723	enddo
724
725	! fabrication de limit.nc
726	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug, &
727	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
728
729
730	call phys_state_var_end
731	if (restart) then
732	print*,'call to restart dyn 1d'
733	Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs, &
734	& u,v,temp,q,omega2)
735
736	print*,'fnday,annee_ref,day_ref,day_ini', &
737	& fnday,annee_ref,day_ref,day_ini
738	day = day_ini
739	day_end = day_ini + nday
740	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
741
742	! Print out the actual date of the beginning of the simulation :
743	call ju2ymds(daytime, an, mois, jour, heure)
744	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
745
746	day = int(daytime)
747	time=daytime-day
748
749	print,'*** intialised fields from restart1dyn *****'
750	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
751	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
752	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
753	! raz for safety
754	do l=1,llm
755	d_q_vert_adv(l,:) = 0.
756	enddo
757	endif
758	!====================== end restart =================================
759	IF (ecrit_slab_oc.eq.1) then
760	open(97,file='div_slab.dat',STATUS='UNKNOWN')
761	elseif (ecrit_slab_oc.eq.0) then
762	open(97,file='div_slab.dat',STATUS='OLD')
763	endif
764	!
765	!=====================================================================
766	IF (CPPKEY_OUTPUTPHYSSCM) THEN
767	CALL iophys_ini(timestep)
768	END IF
769
770	!=====================================================================
771	! START OF THE TEMPORAL LOOP :
772	!=====================================================================
773
774	it_end = nint(fnday*day_step)
775	do while(it.le.it_end)
776
777	if (prt_level.ge.1) then
778	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
779	& it,day,time,it_end,day_step
780	print*,'PAS DE TEMPS ',timestep
781	endif
782	if (it.eq.it_end) lastcall=.True.
783
784	!---------------------------------------------------------------------
785	! Interpolation of forcings in time and onto model levels
786	!---------------------------------------------------------------------
787
788	INCLUDE "1D_interp_cases.h"
789
790	!---------------------------------------------------------------------
791	! Geopotential :
792	!---------------------------------------------------------------------
793	phis(1)=zsurf*RG
794
795	! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
796	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
797
798	do l = 1, llm-1
799	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
800	& (play(l)-play(l+1))/(play(l)+play(l+1))
801	enddo
802
803	!---------------------------------------------------------------------
804	! Vertical advection
805	!---------------------------------------------------------------------
806
807	IF ( forc_w+forc_omega > 0 ) THEN
808
809	IF ( forc_w == 1 ) THEN
810	w_adv=w_mod_cas
811	z_adv=phi/RG
812	ELSE
813	w_adv=omega
814	z_adv=play
815	ENDIF
816
817	! vertical tendencies computed as d X / d t = -W d X / d z
818
819	IF (iflag_1d_vert_adv .EQ. 0) THEN
820
821	! old centered numerical scheme
822	do l=2,llm-1
823	d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
824	d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
825	! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
826	d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l-1))/(z_adv(l+1)-z_adv(l-1))
827	d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1))
828	enddo
829
830
831	ELSE IF (iflag_1d_vert_adv .EQ. 1) THEN
832	! upstream numerical scheme
833
834	do l=2,llm-1
835	IF ( ( ( forc_w .EQ. 1 ) .AND. ( w_adv(l) .GT. 0. ) ) .OR. &
836	( ( forc_w .NE. 1 ) .AND. ( w_adv(l) .LT. 0. ) ) ) THEN
837	d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1))
838	d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1))
839	d_t_vert_adv(l)=-w_adv(l)*(temp(l)-temp(l-1))/(z_adv(l)-z_adv(l-1))
840	d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1))
841	ELSE
842	d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l))
843	d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l))
844	d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l))/(z_adv(l+1)-z_adv(l))
845	d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l))
846	ENDIF
847	enddo
848
849
850	ENDIF ! numerical scheme for vertical advection
851
852
853	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
854	u(1:mxcalc)=u(1:mxcalc) + timestep * d_u_vert_adv(1:mxcalc)
855	v(1:mxcalc)=v(1:mxcalc) + timestep * d_v_vert_adv(1:mxcalc)
856	q(1:mxcalc,:)=q(1:mxcalc,:) + timestep * d_q_vert_adv(1:mxcalc,:)
857	temp(1:mxcalc)=temp(1:mxcalc) + timestep * d_t_vert_adv(1:mxcalc)
858	teta=temp(pzero/play)*rkappa
859	ENDIF
860
861	ENDIF
862
863	!---------------------------------------------------------------------
864	! Listing output for debug prt_level>=1
865	!---------------------------------------------------------------------
866	if (prt_level>=1) then
867	print *,' avant physiq : -------- day time ',day,time
868	write(,) 'firstcall,lastcall,phis', &
869	& firstcall,lastcall,phis
870	end if
871	if (prt_level>=5) then
872	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l', &
873	& 'presniv','plev','play','phi'
874	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l, &
875	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
876	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l', &
877	& 'presniv','u','v','temp','q1','q2','omega2'
878	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l, &
879	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
880	endif
881
882	!---------------------------------------------------------------------
883	! Call physiq :
884	!---------------------------------------------------------------------
885	call physiq(ngrid,llm, &
886	firstcall,lastcall,timestep, &
887	plev,play,phi,phis,presnivs, &
888	u,v, rot, temp,q,omega2, &
889	du_phys,dv_phys,dt_phys,dq,dpsrf)
890	firstcall=.false.
891
892	!---------------------------------------------------------------------
893	! Listing output for debug
894	!---------------------------------------------------------------------
895	if (prt_level>=5) then
896	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l', &
897	& 'presniv','plev','play','phi'
898	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l, &
899	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
900	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l', &
901	& 'presniv','u','v','temp','q1','q2','omega2'
902	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l, &
903	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
904	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l', &
905	& 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
906	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l, &
907	& presnivs(l),86400du_phys(l),86400dv_phys(l), &
908	& 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
909	write(,) 'dpsrf',dpsrf
910	endif
911	!---------------------------------------------------------------------
912	! Add physical tendencies :
913	!---------------------------------------------------------------------
914
915	fcoriolis=2.sin(rpixlat/180.)*romega
916
917	IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
918	fcoriolis, xlat,mxcalc
919
920	!---------------------------------------------------------------------
921	! Geostrophic forcing
922	!---------------------------------------------------------------------
923
924	IF ( forc_geo == 0 ) THEN
925	d_u_age(1:mxcalc)=0.
926	d_v_age(1:mxcalc)=0.
927	ELSE
928	sfdt = sin(0.5fcoriolistimestep)
929	cfdt = cos(0.5fcoriolistimestep)
930
931	d_u_age(1:mxcalc)= -2.sfdt/timestep &
932	& (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - &
933	& cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
934	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
935	!
936	d_v_age(1:mxcalc)= -2.sfdt/timestep &
937	& (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + &
938	& sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
939	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
940	ENDIF
941	!
942	!---------------------------------------------------------------------
943	! Nudging
944	!---------------------------------------------------------------------
945	d_t_nudge(:) = 0.
946	d_u_nudge(:) = 0.
947	d_v_nudge(:) = 0.
948	d_q_nudge(:,:) = 0.
949
950	DO l=1,llm
951
952	IF (nudging_u .LT. 0) THEN
953
954	d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))*invtau_u_nudg_mod_cas(l)
955
956	ELSE
957
958	IF ( play(l) < p_nudging_u .AND. nint(nudging_u) /= 0 ) &
959	& d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))/nudging_u
960
961	ENDIF
962
963
964	IF (nudging_v .LT. 0) THEN
965
966	d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))*invtau_v_nudg_mod_cas(l)
967
968	ELSE
969
970
971	IF ( play(l) < p_nudging_v .AND. nint(nudging_v) /= 0 ) &
972	& d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))/nudging_v
973
974	ENDIF
975
976
977	IF (nudging_t .LT. 0) THEN
978
979	d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))*invtau_temp_nudg_mod_cas(l)
980
981	ELSE
982
983
984	IF ( play(l) < p_nudging_t .AND. nint(nudging_t) /= 0 ) &
985	& d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))/nudging_t
986
987	ENDIF
988
989
990	IF (nudging_qv .LT. 0) THEN
991
992	d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))*invtau_qv_nudg_mod_cas(l)
993
994	ELSE
995
996	IF ( play(l) < p_nudging_qv .AND. nint(nudging_qv) /= 0 ) &
997	& d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv
998
999	ENDIF
1000
1001	ENDDO
1002
1003	!-----------------------------------------------------------
1004	! horizontal forcings (advection) and nudging
1005	!-----------------------------------------------------------
1006
1007	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
1008
1009	u(1:mxcalc)=u(1:mxcalc) + timestep*( &
1010	& du_phys(1:mxcalc) &
1011	& +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) &
1012	& +d_u_nudge(1:mxcalc) )
1013	v(1:mxcalc)=v(1:mxcalc) + timestep*( &
1014	& dv_phys(1:mxcalc) &
1015	& +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) &
1016	& +d_v_nudge(1:mxcalc) )
1017	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( &
1018	& dq(1:mxcalc,:) &
1019	& +d_q_adv(1:mxcalc,:) &
1020	& +d_q_nudge(1:mxcalc,:) )
1021
1022
1023	temp(1:mxcalc)=temp(1:mxcalc)+timestep*( &
1024	& dt_phys(1:mxcalc) &
1025	& +d_t_adv(1:mxcalc) &
1026	& +d_t_nudge(1:mxcalc) &
1027	& +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
1028
1029	!---------------------------------------------------------------------
1030	! Optional outputs
1031	!---------------------------------------------------------------------
1032
1033	IF (CPPKEY_OUTPUTPHYSSCM) THEN
1034	CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
1035	CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv)
1036	CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv)
1037	CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv)
1038	CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1))
1039	CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas)
1040	CALL iophys_ecrit('u',klev,'u','m/s',u)
1041	CALL iophys_ecrit('unud',klev,'unud','m/s',u_nudg_mod_cas)
1042	CALL iophys_ecrit('v',klev,'v','m/s',v)
1043	CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas)
1044	CALL iophys_ecrit('temp',klev,'temp','K',temp)
1045	CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas)
1046	CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge)
1047	CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',100086400d_q_nudge(:,1))
1048	END IF
1049
1050
1051
1052	! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
1053
1054	teta=temp(pzero/play)*rkappa
1055
1056	!---------------------------------------------------------------------
1057	! Nudge soil temperature if requested
1058	!---------------------------------------------------------------------
1059
1060	IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1061	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) &
1062	& -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
1063	ENDIF
1064
1065
1066	END IF ! end if in case tendency should be added
1067
1068	!---------------------------------------------------------------------
1069	! Air temperature :
1070	!---------------------------------------------------------------------
1071	if (lastcall) then
1072	print*,'Pas de temps final ',it
1073	call ju2ymds(daytime, an, mois, jour, heure)
1074	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1075	endif
1076
1077	! incremente day time
1078	daytime = daytime+1./day_step
1079	day = int(daytime+0.1/day_step)
1080	time = time_ini/24.+real(mod(it,day_step))/day_step
1081	it=it+1
1082
1083	enddo
1084
1085	if (ecrit_slab_oc.ne.-1) close(97)
1086
1087	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1088	! ---------------------------------------------------------------------------
1089	call dyn1dredem("restart1dyn.nc", &
1090	& plev,play,phi,phis,presnivs, &
1091	& u,v,temp,q,omega2)
1092
1093	CALL abort_gcm ('lmdz1d ','The End ',0)
1094
1095	END SUBROUTINE scm

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

Download in other formats:

Original Format