Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

old_1D_interp_cases.h

Last change on this file was 5160, checked in by abarral, 5 months ago
Put .h into modules
File size: 39.4 KB

Line
1
2	! $Id: 1D_interp_cases.h 3537 2019-06-19 08:29:16Z fhourdin $
3
4	!---------------------------------------------------------------------
5	! Forcing_LES case: constant dq_dyn
6	!---------------------------------------------------------------------
7	IF (forcing_LES) THEN
8	DO l = 1,llm
9	d_q_adv(l,1) = dq_dyn(l,1)
10	ENDDO
11	endif ! forcing_LES
12	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13	!---------------------------------------------------------------------
14	! Interpolation forcing in time and onto model levels
15	!---------------------------------------------------------------------
16	IF (forcing_GCSSold) THEN
17	CALL get_uvd(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, &
18	& ht_gcssold,hq_gcssold,hw_gcssold, &
19	& hu_gcssold,hv_gcssold, &
20	& hthturb_gcssold,hqturb_gcssold,Ts_gcssold, &
21	& imp_fcg_gcssold,ts_fcg_gcssold, &
22	& Tp_fcg_gcssold,Turb_fcg_gcssold)
23	IF (prt_level.ge.1) THEN
24	PRINT *,' get_uvd -> hqturb_gcssold ',it,hqturb_gcssold
25	endif
26	! large-scale forcing :
27	!!! tsurf = ts_gcssold
28	DO l = 1, llm
29	! u(l) = hu_gcssold(l) ! on prescrit le vent
30	! v(l) = hv_gcssold(l) ! on prescrit le vent
31	! omega(l) = hw_gcssold(l)
32	! rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
33	! omega2(l)=-rho(l)*omega(l)
34	omega(l) = hw_gcssold(l)
35	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
36
37	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
38	d_t_adv(l) = ht_gcssold(l)
39	d_q_adv(l,1) = hq_gcssold(l)
40	dt_cooling(l) = 0.0
41	enddo
42
43	endif ! forcing_GCSSold
44	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
45	!---------------------------------------------------------------------
46	! Interpolation Toga forcing
47	!---------------------------------------------------------------------
48	IF (forcing_toga) THEN
49	IF (prt_level.ge.1) THEN
50	PRINT*, &
51	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_toga=', &
52	& day,day1,(day-day1)86400.,(day-day1)86400/dt_toga
53	endif
54
55	! time interpolation:
56	CALL interp_toga_time(daytime,day1,annee_ref &
57	& ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga &
58	& ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga &
59	& ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga &
60	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
61	& ,ht_prof,vt_prof,hq_prof,vq_prof)
62	! EV: tg instead of ts_cur
63	IF (type_ts_forcing.EQ.1) tg = ts_prof !
64
65	! vertical interpolation:
66	CALL interp_toga_vertical(play,nlev_toga,plev_prof &
67	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
68	& ,ht_prof,vt_prof,hq_prof,vq_prof &
69	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
70	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
71
72	! large-scale forcing :
73	tsurf = ts_prof
74	DO l = 1, llm
75	u(l) = u_mod(l) ! sb: on prescrit le vent
76	v(l) = v_mod(l) ! sb: on prescrit le vent
77	! omega(l) = w_prof(l)
78	! rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
79	! omega2(l)=-rho(l)*omega(l)
80	omega(l) = w_mod(l)
81	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
82
83	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
84	d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
85	d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
86	dt_cooling(l) = 0.0
87	enddo
88
89	endif ! forcing_toga
90	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
91	! Interpolation DICE forcing
92	!---------------------------------------------------------------------
93	IF (forcing_dice) THEN
94	IF (prt_level.ge.1) THEN
95	PRINT,'#### ITAP,day,day1,(day-day1)86400,(day-day1)*86400/dt_dice=',&
96	& day,day1,(day-day1)86400.,(day-day1)86400/dt_dice
97	endif
98
99	! time interpolation:
100	CALL interp_dice_time(daytime,day1,annee_ref &
101	& ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice &
102	& ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice &
103	& ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice &
104	& ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
105	& ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof &
106	& ,ustar_prof,psurf_prof,ug_profd,vg_profd &
107	& ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd &
108	& ,omega_profd)
109	! do l = 1, llm
110	! PRINT *,'llm l omega_profd',llm,l,omega_profd(l)
111	! enddo
112	! EV tg instead of ts_cur
113	IF (type_ts_forcing.EQ.1) tg = tg_prof ! SST used
114
115	! vertical interpolation:
116	CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice &
117	& ,t_dice,qv_dice,u_dice,v_dice,o3_dice &
118	& ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd,omega_profd &
119	& ,t_mod,qv_mod,u_mod,v_mod,o3_mod &
120	& ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
121	! do l = 1, llm
122	! PRINT *,'llm l omega_mod',llm,l,omega_mod(l)
123	! enddo
124
125	! Les forcages DICE sont donnes /jour et non /seconde !
126	ht_mod(:)=ht_mod(:)/86400.
127	hq_mod(:)=hq_mod(:)/86400.
128	hu_mod(:)=hu_mod(:)/86400.
129	hv_mod(:)=hv_mod(:)/86400.
130
131	!calcul de l'advection verticale a partir du omega (repris cas TWPICE, MPL 05082013)
132	!Calcul des gradients verticaux
133	!initialisation
134	d_t_z(:)=0.
135	d_q_z(:)=0.
136	d_u_z(:)=0.
137	d_v_z(:)=0.
138	DO l=2,llm-1
139	d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
140	d_q_z(l)=(q(l+1,1)-q(l-1,1)) /(play(l+1)-play(l-1))
141	d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
142	d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
143	ENDDO
144	d_t_z(1)=d_t_z(2)
145	d_q_z(1)=d_q_z(2)
146	! d_u_z(1)=u(2)/(play(2)-psurf)/5.
147	! d_v_z(1)=v(2)/(play(2)-psurf)/5.
148	d_u_z(1)=0.
149	d_v_z(1)=0.
150	d_t_z(llm)=d_t_z(llm-1)
151	d_q_z(llm)=d_q_z(llm-1)
152	d_u_z(llm)=d_u_z(llm-1)
153	d_v_z(llm)=d_v_z(llm-1)
154
155	!Calcul de l advection verticale:
156	! utiliser omega (Pa/s) et non w (m/s) !! MP 20131108
157	d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
158	d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
159	d_u_dyn_z(:)=omega_mod(:)*d_u_z(:)
160	d_v_dyn_z(:)=omega_mod(:)*d_v_z(:)
161
162	! large-scale forcing :
163	! tsurf = tg_prof MPL 20130925 commente
164	psurf = psurf_prof
165	! For this case, fluxes are imposed
166	fsens=-1*shf_prof
167	flat=-1*lhf_prof
168	ust=ustar_prof
169	tg=tg_prof
170	PRINT *,'ust= ',ust
171	DO l = 1, llm
172	ug(l)= ug_profd
173	vg(l)= vg_profd
174	! omega(l) = w_prof(l)
175	! rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
176	! omega2(l)=-rho(l)*omega(l)
177	! omega(l) = w_mod(l)(-rgrho(l))
178	omega(l) = omega_mod(l)
179	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
180
181	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
182	d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
183	d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
184	d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l)
185	d_v_adv(l) = hv_mod(l)-d_v_dyn_z(l)
186	dt_cooling(l) = 0.0
187	enddo
188
189	endif ! forcing_dice
190	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191	! Interpolation gabls4 forcing
192	!---------------------------------------------------------------------
193	IF (forcing_gabls4 ) THEN
194	IF (prt_level.ge.1) THEN
195	PRINT,'#### ITAP,day,day1,(day-day1)86400,(day-day1)*86400/dt_gabls4=',&
196	& day,day1,(day-day1)86400.,(day-day1)86400/dt_gabls4
197	endif
198
199	! time interpolation:
200	CALL interp_gabls4_time(daytime,day1,annee_ref &
201	& ,year_ini_gabls4,day_ju_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 &
202	& ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 &
203	& ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg)
204	!EV tg instead of ts_cur
205	IF (type_ts_forcing.EQ.1) tg = tg_prof ! SST used
206
207	! vertical interpolation:
208	! on re-utilise le programme interp_dice_vertical: les transformations sur
209	! plev_gabls4,th_gabls4,qv_gabls4,u_gabls4,v_gabls4 ne sont pas prises en compte.
210	! seules celles sur ht_profg,hq_profg,ug_profg,vg_profg sont prises en compte.
211
212	CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 &
213	! & ,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,poub &
214	& ,poub,poub,poub,poub,poub &
215	& ,ht_profg,hq_profg,ug_profg,vg_profg,poub,poub &
216	& ,t_mod,qv_mod,u_mod,v_mod,o3_mod &
217	& ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
218
219	DO l = 1, llm
220	ug(l)= ug_mod(l)
221	vg(l)= vg_mod(l)
222	d_t_adv(l)=ht_mod(l)
223	d_q_adv(l,1)=hq_mod(l)
224	enddo
225
226	endif ! forcing_gabls4
227	!---------------------------------------------------------------------
228
229	!---------------------------------------------------------------------
230	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
231	!---------------------------------------------------------------------
232	! Interpolation forcing TWPice
233	!---------------------------------------------------------------------
234	IF (forcing_twpice) THEN
235	PRINT*, &
236	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_twpi=', &
237	& daytime,day1,(daytime-day1)*86400., &
238	& (daytime-day1)*86400/dt_twpi
239
240	! time interpolation:
241	CALL interp_toga_time(daytime,day1,annee_ref &
242	& ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi &
243	& ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi &
244	& ,ht_twpi,vt_twpi,hq_twpi,vq_twpi &
245	& ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp,u_proftwp &
246	& ,v_proftwp,w_proftwp &
247	& ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
248
249	! vertical interpolation:
250	CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp &
251	& ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp &
252	& ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp &
253	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
254	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
255
256
257	!calcul de l'advection verticale a partir du omega
258	!Calcul des gradients verticaux
259	!initialisation
260	d_t_z(:)=0.
261	d_q_z(:)=0.
262	d_t_dyn_z(:)=0.
263	d_q_dyn_z(:)=0.
264	DO l=2,llm-1
265	d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
266	d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
267	ENDDO
268	d_t_z(1)=d_t_z(2)
269	d_q_z(1)=d_q_z(2)
270	d_t_z(llm)=d_t_z(llm-1)
271	d_q_z(llm)=d_q_z(llm-1)
272
273	!Calcul de l advection verticale
274	d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
275	d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
276
277	!wind nudging above 500m with a 2h time scale
278	DO l=1,llm
279	IF (nudge_wind) THEN
280	! if (phi(l).gt.5000.) THEN
281	IF (phi(l).gt.0.) THEN
282	u(l)=u(l)+timestep(u_mod(l)-u(l))/(2.3600.)
283	v(l)=v(l)+timestep(v_mod(l)-v(l))/(2.3600.)
284	endif
285	else
286	u(l) = u_mod(l)
287	v(l) = v_mod(l)
288	endif
289	enddo
290
291	!CR:nudging of q and theta with a 6h time scale above 15km
292	IF (nudge_thermo) THEN
293	DO l=1,llm
294	zz(l)=phi(l)/9.8
295	IF ((zz(l).le.16000.).AND.(zz(l).gt.15000.)) THEN
296	zfact=(zz(l)-15000.)/1000.
297	q(l,1)=q(l,1)+timestep(q_mod(l)-q(l,1))/(6.3600.)*zfact
298	temp(l)=temp(l)+timestep(t_mod(l)-temp(l))/(6.3600.)*zfact
299	ELSE IF (zz(l).gt.16000.) THEN
300	q(l,1)=q(l,1)+timestep(q_mod(l)-q(l,1))/(6.3600.)
301	temp(l)=temp(l)+timestep(t_mod(l)-temp(l))/(6.3600.)
302	endif
303	enddo
304	endif
305
306	DO l = 1, llm
307	omega(l) = w_mod(l)
308	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
309	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
310	!calcul de l'advection totale
311	IF (cptadvw) THEN
312	d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
313	! PRINT*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
314	d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
315	! PRINT*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
316	else
317	d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
318	d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
319	endif
320	dt_cooling(l) = 0.0
321	enddo
322
323	endif ! forcing_twpice
324
325	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
326	!---------------------------------------------------------------------
327	! Interpolation forcing AMMA
328	!---------------------------------------------------------------------
329
330	IF (forcing_amma) THEN
331	PRINT*, &
332	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_amma=', &
333	& daytime,day1,(daytime-day1)*86400., &
334	& (daytime-day1)*86400/dt_amma
335
336	! time interpolation using TOGA interpolation routine
337	CALL interp_amma_time(daytime,day1,annee_ref &
338	& ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma &
339	& ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma &
340	& ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma &
341	& ,sens_profamma)
342
343	PRINT*,'apres interpolation temporelle AMMA'
344
345	DO k=1,nlev_amma
346	th_profamma(k)=0.
347	q_profamma(k)=0.
348	u_profamma(k)=0.
349	v_profamma(k)=0.
350	vt_profamma(k)=0.
351	vq_profamma(k)=0.
352	enddo
353	! vertical interpolation using TOGA interpolation routine:
354	! WRITE(,)'avant interp vert', t_proftwp
355	CALL interp_toga_vertical(play,nlev_amma,plev_amma &
356	& ,th_profamma,q_profamma,u_profamma,v_profamma &
357	& ,vitw_profamma &
358	& ,ht_profamma,vt_profamma,hq_profamma,vq_profamma &
359	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
360	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
361	WRITE(,) 'Profil initial forcing AMMA interpole'
362
363
364	!calcul de l'advection verticale a partir du omega
365	!Calcul des gradients verticaux
366	!initialisation
367	DO l=1,llm
368	d_t_z(l)=0.
369	d_q_z(l)=0.
370	enddo
371
372	DO l=2,llm-1
373	d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
374	d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
375	ENDDO
376	d_t_z(1)=d_t_z(2)
377	d_q_z(1)=d_q_z(2)
378	d_t_z(llm)=d_t_z(llm-1)
379	d_q_z(llm)=d_q_z(llm-1)
380
381
382	DO l = 1, llm
383	rho(l) = play(l)/(rdtemp(l)(1.+(rv/rd-1.)*q(l,1)))
384	omega(l) = w_mod(l)(-rgrho(l))
385	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
386	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
387	!calcul de l'advection totale
388	! d_t_adv(l) = alphaomega(l)/rcpd+ht_mod(l)-omega(l)d_t_z(l)
389	!attention: on impose dth
390	d_t_adv(l) = alpha*omega(l)/rcpd+ &
391	& ht_mod(l)(play(l)/pzero)rkappa-omega(l)d_t_z(l)
392	! d_t_adv(l) = 0.
393	! PRINT*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
394	d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
395	! d_q_adv(l,1) = 0.
396	! PRINT*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
397
398	dt_cooling(l) = 0.0
399	enddo
400
401
402	! ok_flux_surf=.FALSE.
403	fsens=-1.*sens_profamma
404	flat=-1.*lat_profamma
405
406	endif ! forcing_amma
407
408	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
409	!---------------------------------------------------------------------
410	! Interpolation forcing Rico
411	!---------------------------------------------------------------------
412	IF (forcing_rico) THEN
413	! CALL lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,q,temp,u,v,play)
414	CALL lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
415
416	DO l=1,llm
417	d_t_adv(l) = (dth_rico(l) + dt_dyn(l))
418	d_q_adv(l,1) = (dqh_rico(l) + dq_dyn(l,1))
419	d_q_adv(l,2) = 0.
420	enddo
421	endif ! forcing_rico
422	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
423	!---------------------------------------------------------------------
424	! Interpolation forcing Arm_cu
425	!---------------------------------------------------------------------
426	IF (forcing_armcu) THEN
427	PRINT*, &
428	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_armcu=', &
429	& day,day1,(day-day1)86400.,(day-day1)86400/dt_armcu
430
431	! time interpolation:
432	! ATTENTION, cet appel ne convient pas pour TOGA !!
433	! revoir 1DUTILS.h et les arguments
434	CALL interp_armcu_time(daytime,day1,annee_ref &
435	& ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu &
436	& ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu &
437	& ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof &
438	& ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
439
440	! vertical interpolation:
441	! No vertical interpolation if nlev imposed to 19 or 40
442
443	! For this case, fluxes are imposed
444	fsens=-1*sens_prof
445	flat=-1*flat_prof
446
447	! Advective forcings are given in K or g/kg ... BY HOUR
448	DO l = 1, llm
449	ug(l)= u_mod(l)
450	vg(l)= v_mod(l)
451	IF((phi(l)/RG).LT.1000) THEN
452	d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
453	d_q_adv(l,1) = adv_qt_prof/1000./3600.
454	d_q_adv(l,2) = 0.0
455	! PRINT *,'INF1000: phi dth dq1 dq2',
456	! : phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
457	ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
458	fact=((phi(l)/RG)-1000.)/2000.
459	fact=1-fact
460	d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
461	d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
462	d_q_adv(l,2) = 0.0
463	! PRINT *,'SUP1000: phi fact dth dq1 dq2',
464	! : phi(l)/RG,fact,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
465	ELSE
466	d_t_adv(l) = 0.0
467	d_q_adv(l,1) = 0.0
468	d_q_adv(l,2) = 0.0
469	! PRINT *,'SUP3000: phi dth dq1 dq2',
470	! : phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
471	ENDIF
472	dt_cooling(l) = 0.0
473	! PRINT *,'Interp armcu: phi dth dq1 dq2',
474	! : l,phi(l),d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
475	enddo
476	endif ! forcing_armcu
477	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
478	!---------------------------------------------------------------------
479	! Interpolation forcing in time and onto model levels
480	!---------------------------------------------------------------------
481	IF (forcing_sandu) THEN
482	PRINT*, &
483	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_sandu=', &
484	& day,day1,(day-day1)86400.,(day-day1)86400/dt_sandu
485
486	! time interpolation:
487	! ATTENTION, cet appel ne convient pas pour TOGA !!
488	! revoir 1DUTILS.h et les arguments
489	CALL interp_sandu_time(daytime,day1,annee_ref &
490	& ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu &
491	& ,nlev_sandu &
492	& ,ts_sandu,ts_prof)
493	! EV tg instead of ts_cur
494	IF (type_ts_forcing.EQ.1) tg = ts_prof ! SST used in read_tsurf1d
495
496	! vertical interpolation:
497	CALL interp_sandu_vertical(play,nlev_sandu,plev_profs &
498	& ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs &
499	& ,omega_profs,o3mmr_profs &
500	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
501	& ,omega_mod,o3mmr_mod,mxcalc)
502	!calcul de l'advection verticale
503	!Calcul des gradients verticaux
504	!initialisation
505	d_t_z(:)=0.
506	d_q_z(:)=0.
507	d_t_dyn_z(:)=0.
508	d_q_dyn_z(:)=0.
509	! schema centre
510	! DO l=2,llm-1
511	! d_t_z(l)=(temp(l+1)-temp(l-1))
512	! & /(play(l+1)-play(l-1))
513	! d_q_z(l)=(q(l+1,1)-q(l-1,1))
514	! & /(play(l+1)-play(l-1))
515	! schema amont
516	DO l=2,llm-1
517	d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
518	d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
519	! PRINT *,'l temp2 temp0 play2 play0 omega_mod',
520	! & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
521	ENDDO
522	d_t_z(1)=d_t_z(2)
523	d_q_z(1)=d_q_z(2)
524	d_t_z(llm)=d_t_z(llm-1)
525	d_q_z(llm)=d_q_z(llm-1)
526
527	! calcul de l advection verticale
528	! Confusion w (m/s) et omega (Pa/s) !!
529	d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
530	d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
531	! do l=1,llm
532	! PRINT *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
533	! :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
534	! enddo
535
536
537	! large-scale forcing : pour le cas Sandu ces forcages sont la SST
538	! et une divergence constante -> profil de omega
539	tsurf = ts_prof
540	WRITE(,) 'SST suivante: ',tsurf
541	DO l = 1, llm
542	omega(l) = omega_mod(l)
543	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
544
545	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
546
547	! d_t_adv(l) = 0.0
548	! d_q_adv(l,1) = 0.0
549	!CR:test advection=0
550	!calcul de l'advection verticale
551	d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
552	! PRINT*,'temp adv',l,-d_t_dyn_z(l)
553	d_q_adv(l,1) = -d_q_dyn_z(l)
554	! PRINT*,'q adv',l,-d_q_dyn_z(l)
555	dt_cooling(l) = 0.0
556	enddo
557	endif ! forcing_sandu
558	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
559	!---------------------------------------------------------------------
560	! Interpolation forcing in time and onto model levels
561	!---------------------------------------------------------------------
562	IF (forcing_astex) THEN
563	PRINT*, &
564	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/dt_astex=', &
565	& day,day1,(day-day1)86400.,(day-day1)86400/dt_astex
566
567	! time interpolation:
568	! ATTENTION, cet appel ne convient pas pour TOGA !!
569	! revoir 1DUTILS.h et les arguments
570	CALL interp_astex_time(daytime,day1,annee_ref &
571	& ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex &
572	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
573	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
574	& ,ufa_prof,vfa_prof)
575	! EV tg instead of ts_cur
576	IF (type_ts_forcing.EQ.1) tg = ts_prof ! SST used
577	! vertical interpolation:
578	CALL interp_astex_vertical(play,nlev_astex,plev_profa &
579	& ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa &
580	& ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa &
581	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
582	& ,tke_mod,o3mmr_mod,mxcalc)
583	!calcul de l'advection verticale
584	!Calcul des gradients verticaux
585	!initialisation
586	d_t_z(:)=0.
587	d_q_z(:)=0.
588	d_t_dyn_z(:)=0.
589	d_q_dyn_z(:)=0.
590	! schema centre
591	! DO l=2,llm-1
592	! d_t_z(l)=(temp(l+1)-temp(l-1))
593	! & /(play(l+1)-play(l-1))
594	! d_q_z(l)=(q(l+1,1)-q(l-1,1))
595	! & /(play(l+1)-play(l-1))
596	! schema amont
597	DO l=2,llm-1
598	d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
599	d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
600	! PRINT *,'l temp2 temp0 play2 play0 omega_mod',
601	! & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
602	ENDDO
603	d_t_z(1)=d_t_z(2)
604	d_q_z(1)=d_q_z(2)
605	d_t_z(llm)=d_t_z(llm-1)
606	d_q_z(llm)=d_q_z(llm-1)
607
608	! calcul de l advection verticale
609	! Confusion w (m/s) et omega (Pa/s) !!
610	d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
611	d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
612	! do l=1,llm
613	! PRINT *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
614	! :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
615	! enddo
616
617
618	! large-scale forcing : pour le cas Astex ces forcages sont la SST
619	! la divergence,ug,vg,ufa,vfa
620	tsurf = ts_prof
621	WRITE(,) 'SST suivante: ',tsurf
622	DO l = 1, llm
623	omega(l) = w_mod(l)
624	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
625
626	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
627
628	! d_t_adv(l) = 0.0
629	! d_q_adv(l,1) = 0.0
630	!CR:test advection=0
631	!calcul de l'advection verticale
632	d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
633	! PRINT*,'temp adv',l,-d_t_dyn_z(l)
634	d_q_adv(l,1) = -d_q_dyn_z(l)
635	! PRINT*,'q adv',l,-d_q_dyn_z(l)
636	dt_cooling(l) = 0.0
637	enddo
638	endif ! forcing_astex
639
640	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
641	!---------------------------------------------------------------------
642	! Interpolation forcing standard case
643	!---------------------------------------------------------------------
644	IF (forcing_case) THEN
645	PRINT*,'FORCING CASE forcing_case'
646
647	PRINT*, &
648	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/pdt_cas=', &
649	& daytime,day1,(daytime-day1)*86400., &
650	& (daytime-day1)*86400/pdt_cas
651
652	! time interpolation:
653	CALL interp_case_time(daytime,day1,annee_ref &
654	! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas &
655	& ,nt_cas,nlev_cas &
656	& ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas &
657	& ,vitw_cas,du_cas,hu_cas,vu_cas &
658	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
659	& ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas &
660	& ,uw_cas,vw_cas,q1_cas,q2_cas &
661	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas &
662	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas &
663	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas &
664	& ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas &
665	& ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
666	! EV tg instead of ts_cur
667
668	tg = ts_prof_cas
669	psurf=plev_prof_cas(1)
670
671	! vertical interpolation:
672	CALL interp_case_vertical(play,nlev_cas,plev_prof_cas &
673	& ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &
674	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
675	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
676	& ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &
677	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
678	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
679
680
681	!calcul de l'advection verticale a partir du omega
682	!Calcul des gradients verticaux
683	!initialisation
684	d_t_z(:)=0.
685	d_q_z(:)=0.
686	d_u_z(:)=0.
687	d_v_z(:)=0.
688	d_t_dyn_z(:)=0.
689	d_q_dyn_z(:)=0.
690	d_u_dyn_z(:)=0.
691	d_v_dyn_z(:)=0.
692	DO l=2,llm-1
693	d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
694	d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
695	d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
696	d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
697	ENDDO
698	d_t_z(1)=d_t_z(2)
699	d_q_z(1)=d_q_z(2)
700	d_u_z(1)=d_u_z(2)
701	d_v_z(1)=d_v_z(2)
702	d_t_z(llm)=d_t_z(llm-1)
703	d_q_z(llm)=d_q_z(llm-1)
704	d_u_z(llm)=d_u_z(llm-1)
705	d_v_z(llm)=d_v_z(llm-1)
706
707	!Calcul de l advection verticale
708
709	d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
710
711	d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
712	d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
713	d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
714
715	!wind nudging
716	IF (nudge_u.gt.0.) THEN
717	DO l=1,llm
718	u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
719	enddo
720	else
721	DO l=1,llm
722	u(l) = u_mod_cas(l)
723	enddo
724	endif
725
726	IF (nudge_v.gt.0.) THEN
727	DO l=1,llm
728	v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
729	enddo
730	else
731	DO l=1,llm
732	v(l) = v_mod_cas(l)
733	enddo
734	endif
735
736	IF (nudge_w.gt.0.) THEN
737	DO l=1,llm
738	w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
739	enddo
740	else
741	DO l=1,llm
742	w(l) = w_mod_cas(l)
743	enddo
744	endif
745
746	!nudging of q and temp
747	IF (nudge_t.gt.0.) THEN
748	DO l=1,llm
749	temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
750	enddo
751	endif
752	IF (nudge_q.gt.0.) THEN
753	DO l=1,llm
754	q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
755	enddo
756	endif
757
758	DO l = 1, llm
759	omega(l) = w_mod_cas(l) ! juste car w_mod_cas en Pa/s (MPL 20170310)
760	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
761	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
762
763	!calcul advection
764	IF ((tend_u.EQ.1).AND.(tend_w.EQ.0)) THEN
765	d_u_adv(l)=du_mod_cas(l)
766	ELSE IF ((tend_u.EQ.1).AND.(tend_w.EQ.1)) THEN
767	d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
768	endif
769
770	IF ((tend_v.EQ.1).AND.(tend_w.EQ.0)) THEN
771	d_v_adv(l)=dv_mod_cas(l)
772	ELSE IF ((tend_v.EQ.1).AND.(tend_w.EQ.1)) THEN
773	d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
774	endif
775
776	IF ((tend_t.EQ.1).AND.(tend_w.EQ.0)) THEN
777	! d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
778	d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
779	ELSE IF ((tend_t.EQ.1).AND.(tend_w.EQ.1)) THEN
780	! d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
781	d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
782	endif
783
784	IF ((tend_q.EQ.1).AND.(tend_w.EQ.0)) THEN
785	! d_q_adv(l,1)=dq_mod_cas(l)
786	d_q_adv(l,1)=-1*dq_mod_cas(l)
787	ELSE IF ((tend_q.EQ.1).AND.(tend_w.EQ.1)) THEN
788	! d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
789	d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
790	endif
791
792	IF (tend_rayo.EQ.1) THEN
793	dt_cooling(l) = dtrad_mod_cas(l)
794	! PRINT *,'dt_cooling=',dt_cooling(l)
795	else
796	dt_cooling(l) = 0.0
797	endif
798	enddo
799
800	! Faut-il multiplier par -1 ? (MPL 20160713)
801	IF(ok_flux_surf) THEN
802	fsens=sens_prof_cas
803	flat=lat_prof_cas
804	ENDIF
805
806	IF (ok_prescr_ust) THEN
807	ust=ustar_prof_cas
808	PRINT *,'ust=',ust
809	ENDIF
810	endif ! forcing_case
811
812	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
813	!---------------------------------------------------------------------
814	! Interpolation forcing standard case
815	!---------------------------------------------------------------------
816	IF (forcing_case2 .OR. forcing_SCM) THEN
817	PRINT*,'FORCING CASE forcing_case2'
818	PRINT*, &
819	& '#### ITAP,day,day1,(day-day1)86400,(day-day1)86400/pdt_cas=', &
820	& daytime,day1,(daytime-day1)*86400., &
821	& (daytime-day1)*86400/pdt_cas
822
823	! time interpolation:
824	CALL interp2_case_time(daytime,day1,annee_ref &
825	! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas &
826	& ,nt_cas,nlev_cas &
827	& ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas &
828	& ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas &
829	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
830	& ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas &
831	& ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas &
832
833	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas &
834	& ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas &
835	& ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
836	& ,du_prof_cas,hu_prof_cas,vu_prof_cas &
837	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas &
838	& ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
839	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas &
840	& ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
841	! EV tg instead of ts_cur
842
843	tg = ts_prof_cas
844	! psurf=plev_prof_cas(1)
845	psurf=ps_prof_cas
846
847	! vertical interpolation:
848	CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas &
849	& ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas &
850	& ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas &
851	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
852	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
853	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
854	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas &
855
856	& ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas &
857	& ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas &
858	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
859	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas &
860	& ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
861
862
863	DO l=1,llm
864	teta(l)=temp(l)(100000./play(l))*(rd/rcpd)
865	ENDDO
866	!calcul de l'advection verticale a partir du omega
867	!Calcul des gradients verticaux
868	!initialisation
869	d_t_z(:)=0.
870	d_th_z(:)=0.
871	d_q_z(:)=0.
872	d_u_z(:)=0.
873	d_v_z(:)=0.
874	d_t_dyn_z(:)=0.
875	d_th_dyn_z(:)=0.
876	d_q_dyn_z(:)=0.
877	d_u_dyn_z(:)=0.
878	d_v_dyn_z(:)=0.
879	DO l=2,llm-1
880	d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
881	d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
882	d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
883	d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
884	d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
885	ENDDO
886	d_t_z(1)=d_t_z(2)
887	d_th_z(1)=d_th_z(2)
888	d_q_z(1)=d_q_z(2)
889	d_u_z(1)=d_u_z(2)
890	d_v_z(1)=d_v_z(2)
891	d_t_z(llm)=d_t_z(llm-1)
892	d_th_z(llm)=d_th_z(llm-1)
893	d_q_z(llm)=d_q_z(llm-1)
894	d_u_z(llm)=d_u_z(llm-1)
895	d_v_z(llm)=d_v_z(llm-1)
896
897	!Calcul de l advection verticale
898	! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
899	d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
900	d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
901	d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
902	d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
903	d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
904
905	!geostrophic wind
906	IF (forc_geo.EQ.1) THEN
907	DO l=1,llm
908	ug(l) = ug_mod_cas(l)
909	vg(l) = vg_mod_cas(l)
910	enddo
911	endif
912	!wind nudging
913	IF (nudging_u.gt.0.) THEN
914	DO l=1,llm
915	u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
916	enddo
917	! else
918	! do l=1,llm
919	! u(l) = u_mod_cas(l)
920	! enddo
921	endif
922
923	IF (nudging_v.gt.0.) THEN
924	DO l=1,llm
925	v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
926	enddo
927	! else
928	! do l=1,llm
929	! v(l) = v_mod_cas(l)
930	! enddo
931	endif
932
933	IF (nudging_w.gt.0.) THEN
934	DO l=1,llm
935	w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
936	enddo
937	! else
938	! do l=1,llm
939	! w(l) = w_mod_cas(l)
940	! enddo
941	endif
942
943	!nudging of q and temp
944	IF (nudging_t.gt.0.) THEN
945	DO l=1,llm
946	temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
947	enddo
948	endif
949	IF (nudging_qv.gt.0.) THEN
950	DO l=1,llm
951	q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
952	enddo
953	endif
954
955	DO l = 1, llm
956	! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
957	omega(l) = omega_mod_cas(l)
958	omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
959	alpha = rdtemp(l)(1.+(rv/rd-1.)*q(l,1))/play(l)
960
961	!calcul advections
962	IF ((forc_u.EQ.1).AND.(forc_w.EQ.0)) THEN
963	d_u_adv(l)=du_mod_cas(l)
964	ELSE IF ((forc_u.EQ.1).AND.(forc_w.EQ.1)) THEN
965	d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
966	endif
967
968	IF ((forc_v.EQ.1).AND.(forc_w.EQ.0)) THEN
969	d_v_adv(l)=dv_mod_cas(l)
970	ELSE IF ((forc_v.EQ.1).AND.(forc_w.EQ.1)) THEN
971	d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
972	endif
973
974	! Puisque dth a ete converti en dt, on traite de la meme facon
975	! les flags tadv et thadv
976	IF ((tadv.EQ.1.OR.thadv.EQ.1) .AND. (forc_w.EQ.0)) THEN
977	! d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
978	d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
979	ELSE IF ((tadv.EQ.1.OR.thadv.EQ.1) .AND. (forc_w.EQ.1)) THEN
980	! d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
981	d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
982	endif
983
984	! if ((thadv.EQ.1) .AND. (forc_w.EQ.0)) THEN
985	! d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
986	! d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
987	! ELSE IF ((thadv.EQ.1) .AND. (forc_w.EQ.1)) THEN
988	! d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
989	! d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
990	! endif
991
992	IF ((qadv.EQ.1) .AND. (forc_w.EQ.0)) THEN
993	d_q_adv(l,1)=dq_mod_cas(l)
994	! d_q_adv(l,1)=-1*dq_mod_cas(l)
995	ELSE IF ((qadv.EQ.1) .AND. (forc_w.EQ.1)) THEN
996	d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
997	! d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
998	endif
999
1000	IF (trad.EQ.1) THEN
1001	tend_rayo=1
1002	dt_cooling(l) = dtrad_mod_cas(l)
1003	! PRINT *,'dt_cooling=',dt_cooling(l)
1004	else
1005	dt_cooling(l) = 0.0
1006	endif
1007	enddo
1008
1009	! Faut-il multiplier par -1 ? (MPL 20160713)
1010	IF(ok_flux_surf) THEN
1011	fsens=-1.*sens_prof_cas
1012	flat=-1.*lat_prof_cas
1013	PRINT *,'1D_interp: sens,flat',fsens,flat
1014	ENDIF
1015
1016	IF (ok_prescr_ust) THEN
1017	ust=ustar_prof_cas
1018	PRINT *,'ust=',ust
1019	ENDIF
1020	endif ! forcing_case2
1021	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1022

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

Context Navigation

source: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/old_1D_interp_cases.h

Download in other formats: