Context Navigation

source: LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h @ 2764

Last change on this file since 2764 was 2716, checked in by fhourdin, 8 years ago
Inclusion du cas arm_cu2, avec les nouveaux formats de forçage 1D (Marie-Pierre Lefebvre)
Property svn:keywords set to `Id`
File size: 41.8 KB

Line
1	!
2	! $Id: 1D_read_forc_cases.h 2716 2016-11-28 22:01:20Z lguez $
3	!
4	!----------------------------------------------------------------------
5	! forcing_les = .T. : Impose a constant cooling
6	! forcing_radconv = .T. : Pure radiative-convective equilibrium:
7	!----------------------------------------------------------------------
8
9
10	nq1=0
11	nq2=0
12
13	if (forcing_les .or. forcing_radconv &
14	& .or. forcing_GCSSold .or. forcing_fire) then
15
16	if (forcing_fire) then
17	!----------------------------------------------------------------------
18	!read fire forcings from fire.nc
19	!----------------------------------------------------------------------
20	fich_fire='fire.nc'
21	call read_fire(fich_fire,nlev_fire,nt_fire &
22	& ,height,tttprof,qtprof,uprof,vprof,e12prof &
23	& ,ugprof,vgprof,wfls,dqtdxls &
24	& ,dqtdyls,dqtdtls,thlpcar)
25	write(,) 'Forcing FIRE lu'
26	kmax=120 ! nombre de niveaux dans les profils et forcages
27	else
28	!----------------------------------------------------------------------
29	! Read profiles from files: prof.inp.001 and lscale.inp.001
30	! (repris de readlesfiles)
31	!----------------------------------------------------------------------
32
33	call readprofiles(nlev_max,kmax,nqtot,height, &
34	& tttprof,qtprof,uprof,vprof, &
35	& e12prof,ugprof,vgprof, &
36	& wfls,dqtdxls,dqtdyls,dqtdtls, &
37	& thlpcar,qprof,nq1,nq2)
38	endif
39
40	! compute altitudes of play levels.
41	zlay(1) =zsurf + rdtsurf(psurf-play(1))/(rg*psurf)
42	do l = 2,llm
43	zlay(l) = zlay(l-1)+rdtsurf(psurf-play(1))/(rg*psurf)
44	enddo
45
46	!----------------------------------------------------------------------
47	! Interpolation of the profiles given on the input file to
48	! model levels
49	!----------------------------------------------------------------------
50	zlay(1) = zsurf + rdtsurf(psurf-play(1))/(rg*psurf)
51	do l=1,llm
52	! Above the max altutide of the input file
53
54	if (zlay(l)<height(kmax)) mxcalc=l
55
56	frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
57	ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
58	if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
59	temp(l) = ttt(play(l)/pzero)*rkappa
60	teta(l) = ttt
61	else
62	temp(l) = ttt
63	teta(l) = ttt(pzero/play(l))*rkappa
64	endif
65	print *,' temp,teta ',l,temp(l),teta(l)
66	q(l,1) = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
67	u(l) = uprof(kmax)-frac*( uprof(kmax)- uprof(kmax-1))
68	v(l) = vprof(kmax)-frac*( vprof(kmax)- vprof(kmax-1))
69	ug(l) = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
70	vg(l) = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
71	IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2) &
72	& -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2))
73	omega(l)= wfls(kmax)-frac*( wfls(kmax)- wfls(kmax-1))
74
75	dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
76	dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
77	do k=2,kmax
78	frac = (height(k)-zlay(l))/(height(k)-height(k-1))
79	if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
80	if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
81	ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
82	if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
83	temp(l) = ttt(play(l)/pzero)*rkappa
84	teta(l) = ttt
85	else
86	temp(l) = ttt
87	teta(l) = ttt(pzero/play(l))*rkappa
88	endif
89	print *,' temp,teta ',l,temp(l),teta(l)
90	q(l,1) = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
91	u(l) = uprof(k)-frac*( uprof(k)- uprof(k-1))
92	v(l) = vprof(k)-frac*( vprof(k)- vprof(k-1))
93	ug(l) = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
94	vg(l) = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
95	IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2) &
96	& -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2))
97	omega(l)= wfls(k)-frac*( wfls(k)- wfls(k-1))
98	dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
99	dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
100	elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
101	ttt =tttprof(1)
102	if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
103	temp(l) = ttt(play(l)/pzero)*rkappa
104	teta(l) = ttt
105	else
106	temp(l) = ttt
107	teta(l) = ttt(pzero/play(l))*rkappa
108	endif
109	q(l,1) = qtprof(1)
110	u(l) = uprof(1)
111	v(l) = vprof(1)
112	ug(l) = ugprof(1)
113	vg(l) = vgprof(1)
114	omega(l)= wfls(1)
115	IF (nq2>0) q(l,nq1:nq2)=qprof(1,nq1:nq2)
116	dq_dyn(l,1) =dqtdtls(1)
117	dt_cooling(l)=thlpcar(1)
118	endif
119	enddo
120
121	temp(l)=max(min(temp(l),350.),150.)
122	rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
123	if (l .lt. llm) then
124	zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l))
125	endif
126	omega2(l)=-rho(l)*omega(l)
127	omega(l)= omega(l)(-rgrho(l)) !en Pa/s
128	if (l>1) then
129	if(zlay(l-1)>height(kmax)) then
130	omega(l)=0.0
131	omega2(l)=0.0
132	endif
133	endif
134	if(q(l,1)<0.) q(l,1)=0.0
135	q(l,2) = 0.0
136	enddo
137
138	endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
139	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140	!---------------------------------------------------------------------
141	! Forcing for GCSSold:
142	!---------------------------------------------------------------------
143	if (forcing_GCSSold) then
144	fich_gcssold_ctl = './forcing.ctl'
145	fich_gcssold_dat = './forcing8.dat'
146	call copie(llm,play,psurf,fich_gcssold_ctl)
147	call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, &
148	& ht_gcssold,hq_gcssold,hw_gcssold, &
149	& hu_gcssold,hv_gcssold, &
150	& hthturb_gcssold,hqturb_gcssold,Ts_gcssold, &
151	& imp_fcg_gcssold,ts_fcg_gcssold, &
152	& Tp_fcg_gcssold,Turb_fcg_gcssold)
153	print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
154	endif ! forcing_GCSSold
155	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156	!---------------------------------------------------------------------
157	! Forcing for RICO:
158	!---------------------------------------------------------------------
159	if (forcing_rico) then
160
161	! call writefield_phy('omega', omega,llm+1)
162	fich_rico = 'rico.txt'
163	call read_rico(fich_rico,nlev_rico,ps_rico,play &
164	& ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &
165	& ,dth_rico,dqh_rico)
166	print*, ' on a lu et prepare RICO'
167
168	mxcalc=llm
169	print *, airefi, ' airefi '
170	do l = 1, llm
171	rho(l) = play(l)/(rdt_rico(l)(1.+(rv/rd-1.)*q_rico(l)))
172	temp(l) = t_rico(l)
173	q(l,1) = q_rico(l)
174	q(l,2) = 0.0
175	u(l) = u_rico(l)
176	v(l) = v_rico(l)
177	ug(l)=u_rico(l)
178	vg(l)=v_rico(l)
179	omega(l) = -w_rico(l)*rg
180	omega2(l) = omega(l)/rg*airefi
181	enddo
182	endif
183	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184	!---------------------------------------------------------------------
185	! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) :
186	!---------------------------------------------------------------------
187
188	if (forcing_toga) then
189
190	! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
191	fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
192	CALL read_togacoare(fich_toga,nlev_toga,nt_toga &
193	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
194	& ,ht_toga,vt_toga,hq_toga,vq_toga)
195
196	write(,) 'Forcing TOGA lu'
197
198	! time interpolation for initial conditions:
199	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
200	CALL interp_toga_time(daytime,day1,annee_ref &
201	& ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga &
202	& ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga &
203	& ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga &
204	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
205	& ,ht_prof,vt_prof,hq_prof,vq_prof)
206
207	! vertical interpolation:
208	CALL interp_toga_vertical(play,nlev_toga,plev_prof &
209	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
210	& ,ht_prof,vt_prof,hq_prof,vq_prof &
211	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
212	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
213	write(,) 'Profil initial forcing TOGA interpole'
214
215	! initial and boundary conditions :
216	tsurf = ts_prof
217	write(,) 'SST initiale: ',tsurf
218	do l = 1, llm
219	temp(l) = t_mod(l)
220	q(l,1) = q_mod(l)
221	q(l,2) = 0.0
222	u(l) = u_mod(l)
223	v(l) = v_mod(l)
224	omega(l) = w_mod(l)
225	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
226	!? rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
227	!? omega2(l)=-rho(l)*omega(l)
228	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
229	d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
230	d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
231	d_q_adv(l,2) = 0.0
232	enddo
233
234	endif ! forcing_toga
235	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
236	!---------------------------------------------------------------------
237	! Forcing from TWPICE experiment (Shaocheng et al. 2010) :
238	!---------------------------------------------------------------------
239
240	if (forcing_twpice) then
241	!read TWP-ICE forcings
242	fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
243	call read_twpice(fich_twpice,nlev_twpi,nt_twpi &
244	& ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi &
245	& ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
246
247	write(,) 'Forcing TWP-ICE lu'
248	!Time interpolation for initial conditions using TOGA interpolation routine
249	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1
250	CALL interp_toga_time(daytime,day1,annee_ref &
251	& ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi &
252	& ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi &
253	& ,ht_twpi,vt_twpi,hq_twpi,vq_twpi &
254	& ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp &
255	& ,u_proftwp,v_proftwp,w_proftwp &
256	& ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
257
258	! vertical interpolation using TOGA interpolation routine:
259	! write(,)'avant interp vert', t_proftwp
260	CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp &
261	& ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp &
262	& ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp &
263	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
264	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
265	! write(,) 'Profil initial forcing TWP-ICE interpole',t_mod
266
267	! initial and boundary conditions :
268	! tsurf = ts_proftwp
269	write(,) 'SST initiale: ',tsurf
270	do l = 1, llm
271	temp(l) = t_mod(l)
272	q(l,1) = q_mod(l)
273	q(l,2) = 0.0
274	u(l) = u_mod(l)
275	v(l) = v_mod(l)
276	omega(l) = w_mod(l)
277	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
278
279	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
280	!on applique le forcage total au premier pas de temps
281	!attention: signe different de toga
282	d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
283	d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
284	d_q_adv(l,2) = 0.0
285	enddo
286
287	endif !forcing_twpice
288
289	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290	!---------------------------------------------------------------------
291	! Forcing from AMMA experiment (Couvreux et al. 2010) :
292	!---------------------------------------------------------------------
293
294	if (forcing_amma) then
295
296	call read_1D_cases
297
298	write(,) 'Forcing AMMA lu'
299
300	!champs initiaux:
301	do k=1,nlev_amma
302	th_ammai(k)=th_amma(k)
303	q_ammai(k)=q_amma(k)
304	u_ammai(k)=u_amma(k)
305	v_ammai(k)=v_amma(k)
306	vitw_ammai(k)=vitw_amma(k,12)
307	ht_ammai(k)=ht_amma(k,12)
308	hq_ammai(k)=hq_amma(k,12)
309	vt_ammai(k)=0.
310	vq_ammai(k)=0.
311	enddo
312	omega(:)=0.
313	omega2(:)=0.
314	rho(:)=0.
315	! vertical interpolation using TOGA interpolation routine:
316	! write(,)'avant interp vert', t_proftwp
317	CALL interp_toga_vertical(play,nlev_amma,plev_amma &
318	& ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai &
319	& ,ht_ammai,vt_ammai,hq_ammai,vq_ammai &
320	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
321	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
322	! write(,) 'Profil initial forcing TWP-ICE interpole',t_mod
323
324	! initial and boundary conditions :
325	! tsurf = ts_proftwp
326	write(,) 'SST initiale mxcalc: ',tsurf,mxcalc
327	do l = 1, llm
328	! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
329	! temp(l) = t_mod(l)(play(l)/pzero)*rkappa
330	temp(l) = t_mod(l)
331	q(l,1) = q_mod(l)
332	q(l,2) = 0.0
333	! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
334	u(l) = u_mod(l)
335	v(l) = v_mod(l)
336	rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
337	omega(l) = w_mod(l)(-rgrho(l))
338	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
339
340	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
341	!on applique le forcage total au premier pas de temps
342	!attention: signe different de toga
343	d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
344	!forcage en th
345	! d_th_adv(l) = ht_mod(l)
346	d_q_adv(l,1) = hq_mod(l)
347	d_q_adv(l,2) = 0.0
348	dt_cooling(l)=0.
349	enddo
350	write(,) 'Prof initeforcing AMMA interpole temp39',temp(39)
351
352
353	! ok_flux_surf=.false.
354	fsens=-1.*sens_amma(12)
355	flat=-1.*lat_amma(12)
356
357	endif !forcing_amma
358
359
360	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
361	!---------------------------------------------------------------------
362	! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf)
363	!---------------------------------------------------------------------
364
365	if (forcing_dice) then
366	!read DICE forcings
367	fich_dice='dice_driver.nc'
368	call read_dice(fich_dice,nlev_dice,nt_dice &
369	& ,zz_dice,plev_dice,t_dice,qv_dice,u_dice,v_dice,o3_dice &
370	& ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice&
371	& ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice &
372	& ,hu_dice,hv_dice,w_dice,omega_dice)
373
374	write(,) 'Forcing DICE lu'
375
376	!champs initiaux:
377	do k=1,nlev_dice
378	t_dicei(k)=t_dice(k)
379	qv_dicei(k)=qv_dice(k)
380	u_dicei(k)=u_dice(k)
381	v_dicei(k)=v_dice(k)
382	o3_dicei(k)=o3_dice(k)
383	ht_dicei(k)=ht_dice(k,1)
384	hq_dicei(k)=hq_dice(k,1)
385	hu_dicei(k)=hu_dice(k,1)
386	hv_dicei(k)=hv_dice(k,1)
387	w_dicei(k)=w_dice(k,1)
388	omega_dicei(k)=omega_dice(k,1)
389	enddo
390	omega(:)=0.
391	omega2(:)=0.
392	rho(:)=0.
393	! vertical interpolation using TOGA interpolation routine:
394	! write(,)'avant interp vert', t_proftwp
395	!
396	! CALL interp_dice_time(daytime,day1,annee_ref
397	! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
398	! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
399	! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
400	! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
401	! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
402	! o ,ustar_prof,psurf_prof,ug_profd,vg_profd
403	! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
404	! o ,omega_profd)
405
406	CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice &
407	& ,t_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei &
408	& ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei&
409	& ,t_mod,qv_mod,u_mod,v_mod,o3_mod &
410	& ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
411
412	! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod ?? zero (MPL 20131108)
413	! w_mod(:,:)=0.
414	! omega_mod(:,:)=0.
415
416	! write(,) 'Profil initial forcing DICE interpole',t_mod
417	! Les forcages DICE sont donnes /jour et non /seconde !
418	ht_mod(:)=ht_mod(:)/86400.
419	hq_mod(:)=hq_mod(:)/86400.
420	hu_mod(:)=hu_mod(:)/86400.
421	hv_mod(:)=hv_mod(:)/86400.
422
423	! initial and boundary conditions :
424	write(,) 'SST initiale mxcalc: ',tsurf,mxcalc
425	do l = 1, llm
426	! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
427	! temp(l) = th_mod(l)(play(l)/pzero)*rkappa
428	temp(l) = t_mod(l)
429	q(l,1) = qv_mod(l)
430	q(l,2) = 0.0
431	! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
432	u(l) = u_mod(l)
433	v(l) = v_mod(l)
434	ug(l)=ug_dice(1)
435	vg(l)=vg_dice(1)
436	rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
437	! omega(l) = w_mod(l)(-rgrho(l))
438	omega(l) = omega_mod(l)
439	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
440
441	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
442	!on applique le forcage total au premier pas de temps
443	!attention: signe different de toga
444	d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
445	!forcage en th
446	! d_th_adv(l) = ht_mod(l)
447	d_q_adv(l,1) = hq_mod(l)
448	d_q_adv(l,2) = 0.0
449	dt_cooling(l)=0.
450	enddo
451	write(,) 'Profil initial forcing DICE interpole temp39',temp(39)
452
453
454	! ok_flux_surf=.false.
455	fsens=-1.*shf_dice(1)
456	flat=-1.*lhf_dice(1)
457	! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par
458	! le coefficient de trainee en surface cd*2=ustarvent(k=1)
459	! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
460	! MPL 05082013
461	ust=ustar_dice(1)
462	tg=tg_dice(1)
463	print *,'ust= ',ust
464	IF (tsurf .LE. 0.) THEN
465	tsurf= tg_dice(1)
466	ENDIF
467	psurf= psurf_dice(1)
468	solsw_in = (1.-albedo)/albedo*swup_dice(1)
469	sollw_in = (0.7RSIGMAtemp(1)**4)-lwup_dice(1)
470	PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
471	endif !forcing_dice
472
473	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
474	!---------------------------------------------------------------------
475	! Forcing from GABLS4 experiment
476	!---------------------------------------------------------------------
477
478	!!!! Si la temperature de surface n'est pas impos??e:
479
480	if (forcing_gabls4) then
481	!read GABLS4 forcings
482
483	fich_gabls4='gabls4_driver.nc'
484
485
486	call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 &
487	& ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4)
488
489	write(,) 'Forcing GABLS4 lu'
490
491	!champs initiaux:
492	do k=1,nlev_gabls4
493	t_gabi(k)=t_gabls4(k)
494	qv_gabi(k)=qv_gabls4(k)
495	u_gabi(k)=u_gabls4(k)
496	v_gabi(k)=v_gabls4(k)
497	poub(k)=0.
498	ht_gabi(k)=ht_gabls4(k,1)
499	hq_gabi(k)=hq_gabls4(k,1)
500	ug_gabi(k)=ug_gabls4(k,1)
501	vg_gabi(k)=vg_gabls4(k,1)
502	enddo
503
504	omega(:)=0.
505	omega2(:)=0.
506	rho(:)=0.
507	! vertical interpolation using TOGA interpolation routine:
508	! write(,)'avant interp vert', t_proftwp
509	!
510	! CALL interp_dice_time(daytime,day1,annee_ref
511	! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
512	! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
513	! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
514	! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
515	! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
516	! o ,ustar_prof,psurf_prof,ug_profd,vg_profd
517	! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
518	! o ,omega_profd)
519
520	CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 &
521	& ,t_gabi,qv_gabi,u_gabi,v_gabi,poub &
522	& ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub &
523	& ,t_mod,qv_mod,u_mod,v_mod,o3_mod &
524	& ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
525
526	! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024
527	! ht_mod(:)=ht_mod(:)/86400.
528	! hq_mod(:)=hq_mod(:)/86400.
529
530	! initial and boundary conditions :
531	write(,) 'SST initiale mxcalc: ',tsurf,mxcalc
532	do l = 1, llm
533	! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
534	! temp(l) = th_mod(l)(play(l)/pzero)*rkappa
535	temp(l) = t_mod(l)
536	q(l,1) = qv_mod(l)
537	q(l,2) = 0.0
538	! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
539	u(l) = u_mod(l)
540	v(l) = v_mod(l)
541	ug(l)=ug_mod(l)
542	vg(l)=vg_mod(l)
543
544	!
545	! tg=tsurf
546	!
547
548	print ,'**** tsurf=',tsurf
549	rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
550	! omega(l) = w_mod(l)(-rgrho(l))
551	omega(l) = omega_mod(l)
552	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
553
554
555
556	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
557	!on applique le forcage total au premier pas de temps
558	!attention: signe different de toga
559	! d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
560	!forcage en th
561	d_th_adv(l) = ht_mod(l)
562	d_q_adv(l,1) = hq_mod(l)
563	d_q_adv(l,2) = 0.0
564	dt_cooling(l)=0.
565	enddo
566
567	!--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024---------------
568	! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par
569	! le coefficient de trainee en surface cd*2=ustarvent(k=1)
570	! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
571	! MPL 05082013
572	! ust=ustar_dice(1)
573	! tg=tg_dice(1)
574	! print *,'ust= ',ust
575	! IF (tsurf .LE. 0.) THEN
576	! tsurf= tg_dice(1)
577	! ENDIF
578	! psurf= psurf_dice(1)
579	! solsw_in = (1.-albedo)/albedo*swup_dice(1)
580	! sollw_in = (0.7RSIGMAtemp(1)**4)-lwup_dice(1)
581	! PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
582	!--------------------------------------------------------------------------------------
583	endif !forcing_gabls4
584
585
586
587	! Forcing from Arm_Cu case
588	! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
589	! large scale advective forcing,radiative forcing
590	! and advective tendency of theta and qt to be applied
591	!---------------------------------------------------------------------
592
593	if (forcing_armcu) then
594	! read armcu forcing :
595	write(,) 'Avant lecture Forcing Arm_Cu'
596	fich_armcu = './ifa_armcu.txt'
597	CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu, &
598	& sens_armcu,flat_armcu,adv_theta_armcu, &
599	& rad_theta_armcu,adv_qt_armcu)
600	write(,) 'Forcing Arm_Cu lu'
601
602	!----------------------------------------------------------------------
603	! Read profiles from file: prof.inp.19 or prof.inp.40
604	! For this case, profiles are given for two vertical resolution
605	! 19 or 40 levels
606	!
607	! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html
608	! Note that the initial profiles contain no liquid water!
609	! (so potential temperature can be interpreted as liquid water
610	! potential temperature and water vapor as total water)
611	! profiles are given at full levels
612	!----------------------------------------------------------------------
613
614	call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod, &
615	& v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
616
617	! time interpolation for initial conditions:
618	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
619
620	print *,'Avant interp_armcu_time'
621	print *,'daytime=',daytime
622	print *,'day1=',day1
623	print *,'annee_ref=',annee_ref
624	print *,'year_ini_armcu=',year_ini_armcu
625	print *,'day_ju_ini_armcu=',day_ju_ini_armcu
626	print *,'nt_armcu=',nt_armcu
627	print *,'dt_armcu=',dt_armcu
628	print *,'nlev_armcu=',nlev_armcu
629	CALL interp_armcu_time(daytime,day1,annee_ref &
630	& ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu &
631	& ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu &
632	& ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof &
633	& ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
634	write(,) 'Forcages interpoles dans temps'
635
636	! No vertical interpolation if nlev imposed to 19 or 40
637	! The vertical grid stops at 4000m # 600hPa
638	mxcalc=llm
639
640	! initial and boundary conditions :
641	! tsurf = ts_prof
642	! tsurf read in lmdz1d.def
643	write(,) 'Tsurf initiale: ',tsurf
644	do l = 1, llm
645	play(l)=play_mod(l)*100.
646	presnivs(l)=play(l)
647	zlay(l)=height(l)
648	temp(l) = t_mod(l)
649	teta(l)=theta_mod(l)
650	q(l,1) = qv_mod(l)/1000.
651	! No liquid water in the initial profil
652	q(l,2) = 0.
653	u(l) = u_mod(l)
654	ug(l)= u_mod(l)
655	v(l) = v_mod(l)
656	vg(l)= v_mod(l)
657	! Advective forcings are given in K or g/kg ... per HOUR
658	! IF(height(l).LT.1000) THEN
659	! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
660	! d_q_adv(l,1) = adv_qt_prof/1000./3600.
661	! d_q_adv(l,2) = 0.0
662	! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
663	! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
664	! : (1-(height(l)-1000.)/2000.)
665	! d_th_adv(l) = d_th_adv(l)/3600.
666	! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
667	! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
668	! d_q_adv(l,2) = 0.0
669	! ELSE
670	! d_th_adv(l) = 0.0
671	! d_q_adv(l,1) = 0.0
672	! d_q_adv(l,2) = 0.0
673	! ENDIF
674	enddo
675	! plev at half levels is given in proh.inp.19 or proh.inp.40 files
676	plev(1)= ap(llm+1)+bp(llm+1)*psurf
677	do l = 1, llm
678	plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
679	print *,'Read_forc: l height play plev zlay temp', &
680	& l,height(l),play(l),plev(l),zlay(l),temp(l)
681	enddo
682	! For this case, fluxes are imposed
683	fsens=-1*sens_prof
684	flat=-1*flat_prof
685
686	endif ! forcing_armcu
687	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
688	!---------------------------------------------------------------------
689	! Forcing from transition case of Irina Sandu
690	!---------------------------------------------------------------------
691
692	if (forcing_sandu) then
693	write(,) 'Avant lecture Forcing SANDU'
694
695	! read sanduref forcing :
696	fich_sandu = './ifa_sanduref.txt'
697	CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
698
699	write(,) 'Forcing SANDU lu'
700
701	!----------------------------------------------------------------------
702	! Read profiles from file: prof.inp.001
703	!----------------------------------------------------------------------
704
705	call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs, &
706	& thl_profs,q_profs,u_profs,v_profs, &
707	& w_profs,omega_profs,o3mmr_profs)
708
709	! time interpolation for initial conditions:
710	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
711	! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
712	! revoir 1DUTILS.h et les arguments
713
714	print *,'Avant interp_sandu_time'
715	print *,'daytime=',daytime
716	print *,'day1=',day1
717	print *,'annee_ref=',annee_ref
718	print *,'year_ini_sandu=',year_ini_sandu
719	print *,'day_ju_ini_sandu=',day_ju_ini_sandu
720	print *,'nt_sandu=',nt_sandu
721	print *,'dt_sandu=',dt_sandu
722	print *,'nlev_sandu=',nlev_sandu
723	CALL interp_sandu_time(daytime,day1,annee_ref &
724	& ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu &
725	& ,nlev_sandu &
726	& ,ts_sandu,ts_prof)
727
728	! vertical interpolation:
729	print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
730	CALL interp_sandu_vertical(play,nlev_sandu,plev_profs &
731	& ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs &
732	& ,omega_profs,o3mmr_profs &
733	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
734	& ,omega_mod,o3mmr_mod,mxcalc)
735	write(,) 'Profil initial forcing SANDU interpole'
736
737	! initial and boundary conditions :
738	tsurf = ts_prof
739	write(,) 'SST initiale: ',tsurf
740	do l = 1, llm
741	temp(l) = t_mod(l)
742	tetal(l)=thl_mod(l)
743	q(l,1) = q_mod(l)
744	q(l,2) = 0.0
745	u(l) = u_mod(l)
746	v(l) = v_mod(l)
747	w(l) = w_mod(l)
748	omega(l) = omega_mod(l)
749	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
750	!? rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
751	!? omega2(l)=-rho(l)*omega(l)
752	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
753	! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
754	! d_q_adv(l,1) = vq_mod(l)
755	d_th_adv(l) = alpha*omega(l)/rcpd
756	d_q_adv(l,1) = 0.0
757	d_q_adv(l,2) = 0.0
758	enddo
759
760	endif ! forcing_sandu
761	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
762	!---------------------------------------------------------------------
763	! Forcing from Astex case
764	!---------------------------------------------------------------------
765
766	if (forcing_astex) then
767	write(,) 'Avant lecture Forcing Astex'
768
769	! read astex forcing :
770	fich_astex = './ifa_astex.txt'
771	CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex, &
772	& ug_astex,vg_astex,ufa_astex,vfa_astex)
773
774	write(,) 'Forcing Astex lu'
775
776	!----------------------------------------------------------------------
777	! Read profiles from file: prof.inp.001
778	!----------------------------------------------------------------------
779
780	call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa, &
781	& thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa, &
782	& w_profa,tke_profa,o3mmr_profa)
783
784	! time interpolation for initial conditions:
785	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
786	! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
787	! revoir 1DUTILS.h et les arguments
788
789	print *,'Avant interp_astex_time'
790	print *,'daytime=',daytime
791	print *,'day1=',day1
792	print *,'annee_ref=',annee_ref
793	print *,'year_ini_astex=',year_ini_astex
794	print *,'day_ju_ini_astex=',day_ju_ini_astex
795	print *,'nt_astex=',nt_astex
796	print *,'dt_astex=',dt_astex
797	print *,'nlev_astex=',nlev_astex
798	CALL interp_astex_time(daytime,day1,annee_ref &
799	& ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex &
800	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
801	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
802	& ,ufa_prof,vfa_prof)
803
804	! vertical interpolation:
805	print *,'Avant interp_vertical: nlev_astex=',nlev_astex
806	CALL interp_astex_vertical(play,nlev_astex,plev_profa &
807	& ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa &
808	& ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa &
809	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
810	& ,tke_mod,o3mmr_mod,mxcalc)
811	write(,) 'Profil initial forcing Astex interpole'
812
813	! initial and boundary conditions :
814	tsurf = ts_prof
815	write(,) 'SST initiale: ',tsurf
816	do l = 1, llm
817	temp(l) = t_mod(l)
818	tetal(l)=thl_mod(l)
819	q(l,1) = qv_mod(l)
820	q(l,2) = ql_mod(l)
821	u(l) = u_mod(l)
822	v(l) = v_mod(l)
823	w(l) = w_mod(l)
824	omega(l) = w_mod(l)
825	! omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
826	! rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
827	! omega2(l)=-rho(l)*omega(l)
828	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
829	! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
830	! d_q_adv(l,1) = vq_mod(l)
831	d_th_adv(l) = alpha*omega(l)/rcpd
832	d_q_adv(l,1) = 0.0
833	d_q_adv(l,2) = 0.0
834	enddo
835
836	endif ! forcing_astex
837	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
838	!---------------------------------------------------------------------
839	! Forcing from standard case :
840	!---------------------------------------------------------------------
841
842	if (forcing_case) then
843
844	write(,),'avant call read_1D_cas'
845	call read_1D_cas
846	write(,) 'Forcing read'
847
848	!Time interpolation for initial conditions using TOGA interpolation routine
849	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1
850	CALL interp_case_time(day,day1,annee_ref &
851	! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas &
852	& ,nt_cas,nlev_cas &
853	& ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas &
854	& ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas &
855	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
856	& ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas &
857	& ,uw_cas,vw_cas,q1_cas,q2_cas &
858	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas &
859	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas &
860	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas &
861	& ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas,ustar_prof_cas &
862	& ,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
863
864	! vertical interpolation using TOGA interpolation routine:
865	! write(,)'avant interp vert', t_prof
866	CALL interp_case_vertical(play,nlev_cas,plev_prof_cas &
867	& ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &
868	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
869	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
870	& ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &
871	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
872	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
873	! write(,) 'Profil initial forcing case interpole',t_mod
874
875	! initial and boundary conditions :
876	! tsurf = ts_prof_cas
877	ts_cur = ts_prof_cas
878	psurf=plev_prof_cas(1)
879	write(,) 'SST initiale: ',tsurf
880	do l = 1, llm
881	temp(l) = t_mod_cas(l)
882	q(l,1) = q_mod_cas(l)
883	q(l,2) = 0.0
884	u(l) = u_mod_cas(l)
885	v(l) = v_mod_cas(l)
886	omega(l) = w_mod_cas(l)
887	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
888
889	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
890	!on applique le forcage total au premier pas de temps
891	!attention: signe different de toga
892	d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
893	d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
894	d_q_adv(l,2) = 0.0
895	d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
896	d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
897	enddo
898
899	! In case fluxes are imposed
900	IF (ok_flux_surf) THEN
901	fsens=sens_prof_cas
902	flat=lat_prof_cas
903	ENDIF
904	IF (ok_prescr_ust) THEN
905	ust=ustar_prof_cas
906	print *,'ust=',ust
907	ENDIF
908
909	endif !forcing_case
910	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
911	!---------------------------------------------------------------------
912	! Forcing from standard case :
913	!---------------------------------------------------------------------
914
915	if (forcing_case2) then
916
917	write(,),'avant call read2_1D_cas'
918	call read2_1D_cas
919	write(,) 'Forcing read'
920
921	!Time interpolation for initial conditions using interpolation routine
922	write(,) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1
923	CALL interp2_case_time(daytime,day1,annee_ref &
924	! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas &
925	& ,nt_cas,nlev_cas &
926	& ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas &
927	& ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas &
928	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
929	& ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas &
930	& ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas &
931	!
932	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas &
933	& ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas &
934	& ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
935	& ,du_prof_cas,hu_prof_cas,vu_prof_cas &
936	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas &
937	& ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
938	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas &
939	& ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
940
941	do l = 1, nlev_cas
942	print *,'apres 1ere interp: plev_cas, plev_prof_cas=',l,plev_cas(l,1),plev_prof_cas(l)
943	enddo
944
945	! vertical interpolation using interpolation routine:
946	! write(,)'avant interp vert', t_prof
947	CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas &
948	& ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas &
949	& ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas &
950	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
951	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
952	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
953	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas &
954	!
955	& ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas &
956	& ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas &
957	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
958	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas &
959	& ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
960
961	! write(,) 'Profil initial forcing case interpole',t_mod
962
963	! initial and boundary conditions :
964	! tsurf = ts_prof_cas
965	ts_cur = ts_prof_cas
966	psurf=plev_prof_cas(1)
967	write(,) 'SST initiale: ',tsurf
968	do l = 1, llm
969	temp(l) = t_mod_cas(l)
970	q(l,1) = qv_mod_cas(l)
971	q(l,2) = ql_mod_cas(l)
972	u(l) = u_mod_cas(l)
973	ug(l)= u_mod_cas(l)
974	v(l) = v_mod_cas(l)
975	vg(l)= v_mod_cas(l)
976	omega(l) = w_mod_cas(l)
977	omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
978
979	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
980	!on applique le forcage total au premier pas de temps
981	!attention: signe different de toga
982	d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
983	d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
984	! d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
985	d_q_adv(l,1) = dq_mod_cas(l)
986	d_q_adv(l,2) = 0.0
987	! d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
988	d_u_adv(l) = du_mod_cas(l)
989	! d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
990	d_u_adv(l) = dv_mod_cas(l)
991	enddo
992
993	! Faut-il multiplier par -1 ? (MPL 20160713)
994	IF (ok_flux_surf) THEN
995	fsens=-1.*sens_prof_cas
996	flat=-1.*lat_prof_cas
997	ENDIF
998	!
999	IF (ok_prescr_ust) THEN
1000	ust=ustar_prof_cas
1001	print *,'ust=',ust
1002	ENDIF
1003
1004	endif !forcing_case2
1005	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006

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

Download in other formats: