Context Navigation

mod_1D_cases_read2.F90 @ 2740

Last change on this file since 2740 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)
File size: 49.1 KB

Line
1	!
2	! $Id: mod_1D_cases_read.F90 2373 2015-10-13 17:28:01Z jyg $
3	!
4	MODULE mod_1D_cases_read2
5
6	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7	!Declarations specifiques au cas standard
8	character*80 :: fich_cas
9	! Discr?tisation
10	integer nlev_cas, nt_cas
11
12
13	!profils environnementaux
14	real, allocatable:: plev_cas(:,:),plevh_cas(:)
15	real, allocatable:: ap_cas(:),bp_cas(:)
16
17	real, allocatable:: z_cas(:,:),zh_cas(:)
18	real, allocatable:: t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
19	real, allocatable:: th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
20	real, allocatable:: u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:)
21
22	!forcing
23	real, allocatable:: ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
24	real, allocatable:: hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
25	real, allocatable:: hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
26	real, allocatable:: hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
27	real, allocatable:: hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
28	real, allocatable:: hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
29	real, allocatable:: ug_cas(:,:),vg_cas(:,:)
30	real, allocatable:: lat_cas(:),sens_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
31	real, allocatable:: uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tke_cas(:)
32
33	!champs interpoles
34	real, allocatable:: plev_prof_cas(:)
35	real, allocatable:: t_prof_cas(:)
36	real, allocatable:: theta_prof_cas(:)
37	real, allocatable:: thl_prof_cas(:)
38	real, allocatable:: thv_prof_cas(:)
39	real, allocatable:: q_prof_cas(:)
40	real, allocatable:: qv_prof_cas(:)
41	real, allocatable:: ql_prof_cas(:)
42	real, allocatable:: qi_prof_cas(:)
43	real, allocatable:: rh_prof_cas(:)
44	real, allocatable:: rv_prof_cas(:)
45	real, allocatable:: u_prof_cas(:)
46	real, allocatable:: v_prof_cas(:)
47	real, allocatable:: vitw_prof_cas(:)
48	real, allocatable:: omega_prof_cas(:)
49	real, allocatable:: ug_prof_cas(:)
50	real, allocatable:: vg_prof_cas(:)
51	real, allocatable:: ht_prof_cas(:)
52	real, allocatable:: hth_prof_cas(:)
53	real, allocatable:: hq_prof_cas(:)
54	real, allocatable:: vt_prof_cas(:)
55	real, allocatable:: vth_prof_cas(:)
56	real, allocatable:: vq_prof_cas(:)
57	real, allocatable:: dt_prof_cas(:)
58	real, allocatable:: dth_prof_cas(:)
59	real, allocatable:: dtrad_prof_cas(:)
60	real, allocatable:: dq_prof_cas(:)
61	real, allocatable:: hu_prof_cas(:)
62	real, allocatable:: hv_prof_cas(:)
63	real, allocatable:: vu_prof_cas(:)
64	real, allocatable:: vv_prof_cas(:)
65	real, allocatable:: du_prof_cas(:)
66	real, allocatable:: dv_prof_cas(:)
67	real, allocatable:: uw_prof_cas(:)
68	real, allocatable:: vw_prof_cas(:)
69	real, allocatable:: q1_prof_cas(:)
70	real, allocatable:: q2_prof_cas(:)
71
72
73	real lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tke_prof_cas
74	real o3_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
75
76
77
78	CONTAINS
79
80	SUBROUTINE read_1D_cas
81	implicit none
82
83	#include "netcdf.inc"
84
85	INTEGER nid,rid,ierr
86	INTEGER ii,jj
87
88	fich_cas='setup/cas.nc'
89	print*,'fich_cas ',fich_cas
90	ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
91	print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
92	if (ierr.NE.NF_NOERR) then
93	write(,) 'ERROR: GROS Pb opening forcings nc file '
94	write(,) NF_STRERROR(ierr)
95	stop ""
96	endif
97	!.......................................................................
98	ierr=NF_INQ_DIMID(nid,'lat',rid)
99	IF (ierr.NE.NF_NOERR) THEN
100	print*, 'Oh probleme lecture dimension lat'
101	ENDIF
102	ierr=NF_INQ_DIMLEN(nid,rid,ii)
103	print*,'OK1 nid,rid,lat',nid,rid,ii
104	!.......................................................................
105	ierr=NF_INQ_DIMID(nid,'lon',rid)
106	IF (ierr.NE.NF_NOERR) THEN
107	print*, 'Oh probleme lecture dimension lon'
108	ENDIF
109	ierr=NF_INQ_DIMLEN(nid,rid,jj)
110	print*,'OK2 nid,rid,lat',nid,rid,jj
111	!.......................................................................
112	ierr=NF_INQ_DIMID(nid,'lev',rid)
113	IF (ierr.NE.NF_NOERR) THEN
114	print*, 'Oh probleme lecture dimension zz'
115	ENDIF
116	ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
117	print*,'OK3 nid,rid,nlev_cas',nid,rid,nlev_cas
118	!.......................................................................
119	ierr=NF_INQ_DIMID(nid,'time',rid)
120	print*,'nid,rid',nid,rid
121	nt_cas=0
122	IF (ierr.NE.NF_NOERR) THEN
123	stop 'probleme lecture dimension sens'
124	ENDIF
125	ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
126	print*,'OK4 nid,rid,nt_cas',nid,rid,nt_cas
127
128	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129	!profils moyens:
130	allocate(plev_cas(nlev_cas,nt_cas))
131	allocate(z_cas(nlev_cas,nt_cas))
132	allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
133	allocate(th_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
134	allocate(u_cas(nlev_cas,nt_cas))
135	allocate(v_cas(nlev_cas,nt_cas))
136
137	!forcing
138	allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
139	allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
140	allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
141	allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
142	allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
143	allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
144	allocate(vitw_cas(nlev_cas,nt_cas))
145	allocate(ug_cas(nlev_cas,nt_cas))
146	allocate(vg_cas(nlev_cas,nt_cas))
147	allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas))
148	allocate(uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas),q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas))
149
150
151	!champs interpoles
152	allocate(plev_prof_cas(nlev_cas))
153	allocate(t_prof_cas(nlev_cas))
154	allocate(q_prof_cas(nlev_cas))
155	allocate(u_prof_cas(nlev_cas))
156	allocate(v_prof_cas(nlev_cas))
157
158	allocate(vitw_prof_cas(nlev_cas))
159	allocate(ug_prof_cas(nlev_cas))
160	allocate(vg_prof_cas(nlev_cas))
161	allocate(ht_prof_cas(nlev_cas))
162	allocate(hq_prof_cas(nlev_cas))
163	allocate(hu_prof_cas(nlev_cas))
164	allocate(hv_prof_cas(nlev_cas))
165	allocate(vt_prof_cas(nlev_cas))
166	allocate(vq_prof_cas(nlev_cas))
167	allocate(vu_prof_cas(nlev_cas))
168	allocate(vv_prof_cas(nlev_cas))
169	allocate(dt_prof_cas(nlev_cas))
170	allocate(dtrad_prof_cas(nlev_cas))
171	allocate(dq_prof_cas(nlev_cas))
172	allocate(du_prof_cas(nlev_cas))
173	allocate(dv_prof_cas(nlev_cas))
174	allocate(uw_prof_cas(nlev_cas))
175	allocate(vw_prof_cas(nlev_cas))
176	allocate(q1_prof_cas(nlev_cas))
177	allocate(q2_prof_cas(nlev_cas))
178
179	print*,'Allocations OK'
180	call read_cas(nid,nlev_cas,nt_cas &
181	& ,z_cas,plev_cas,t_cas,q_cas,rh_cas,th_cas,rv_cas,u_cas,v_cas &
182	& ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas,dv_cas,hv_cas,vv_cas &
183	& ,dt_cas,dtrad_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas &
184	& ,dth_cas,hth_cas,vth_cas,dr_cas,hr_cas,vr_cas,sens_cas,lat_cas,ts_cas&
185	& ,ustar_cas,uw_cas,vw_cas,q1_cas,q2_cas)
186	print*,'Read cas OK'
187
188
189	END SUBROUTINE read_1D_cas
190	!**********************************************************************************************
191	SUBROUTINE read2_1D_cas
192	implicit none
193
194	#include "netcdf.inc"
195
196	INTEGER nid,rid,ierr
197	INTEGER ii,jj
198
199	fich_cas='setup/cas.nc'
200	print*,'fich_cas ',fich_cas
201	ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
202	print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
203	if (ierr.NE.NF_NOERR) then
204	write(,) 'ERROR: GROS Pb opening forcings nc file '
205	write(,) NF_STRERROR(ierr)
206	stop ""
207	endif
208	!.......................................................................
209	ierr=NF_INQ_DIMID(nid,'lat',rid)
210	IF (ierr.NE.NF_NOERR) THEN
211	print*, 'Oh probleme lecture dimension lat'
212	ENDIF
213	ierr=NF_INQ_DIMLEN(nid,rid,ii)
214	print*,'OK1 read2: nid,rid,lat',nid,rid,ii
215	!.......................................................................
216	ierr=NF_INQ_DIMID(nid,'lon',rid)
217	IF (ierr.NE.NF_NOERR) THEN
218	print*, 'Oh probleme lecture dimension lon'
219	ENDIF
220	ierr=NF_INQ_DIMLEN(nid,rid,jj)
221	print*,'OK2 read2: nid,rid,lat',nid,rid,jj
222	!.......................................................................
223	ierr=NF_INQ_DIMID(nid,'nlev',rid)
224	IF (ierr.NE.NF_NOERR) THEN
225	print*, 'Oh probleme lecture dimension nlev'
226	ENDIF
227	ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
228	print*,'OK3 read2: nid,rid,nlev_cas',nid,rid,nlev_cas
229	!.......................................................................
230	ierr=NF_INQ_DIMID(nid,'time',rid)
231	nt_cas=0
232	IF (ierr.NE.NF_NOERR) THEN
233	stop 'Oh probleme lecture dimension time'
234	ENDIF
235	ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
236	print*,'OK4 read2: nid,rid,nt_cas',nid,rid,nt_cas
237
238	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
239	!profils moyens:
240	allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))
241	allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
242	allocate(ap_cas(nlev_cas+1),bp_cas(nt_cas+1))
243	allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
244	allocate(th_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
245	allocate(u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas),vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas))
246
247	!forcing
248	allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
249	allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
250	allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
251	allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
252	allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
253	allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
254	allocate(ug_cas(nlev_cas,nt_cas))
255	allocate(vg_cas(nlev_cas,nt_cas))
256	allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tke_cas(nt_cas))
257	allocate(uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas),q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas))
258
259
260
261	!champs interpoles
262	allocate(plev_prof_cas(nlev_cas))
263	allocate(t_prof_cas(nlev_cas))
264	allocate(theta_prof_cas(nlev_cas))
265	allocate(thl_prof_cas(nlev_cas))
266	allocate(thv_prof_cas(nlev_cas))
267	allocate(q_prof_cas(nlev_cas))
268	allocate(qv_prof_cas(nlev_cas))
269	allocate(ql_prof_cas(nlev_cas))
270	allocate(qi_prof_cas(nlev_cas))
271	allocate(rh_prof_cas(nlev_cas))
272	allocate(rv_prof_cas(nlev_cas))
273	allocate(u_prof_cas(nlev_cas))
274	allocate(v_prof_cas(nlev_cas))
275	allocate(vitw_prof_cas(nlev_cas))
276	allocate(omega_prof_cas(nlev_cas))
277	allocate(ug_prof_cas(nlev_cas))
278	allocate(vg_prof_cas(nlev_cas))
279	allocate(ht_prof_cas(nlev_cas))
280	allocate(hth_prof_cas(nlev_cas))
281	allocate(hq_prof_cas(nlev_cas))
282	allocate(hu_prof_cas(nlev_cas))
283	allocate(hv_prof_cas(nlev_cas))
284	allocate(vt_prof_cas(nlev_cas))
285	allocate(vth_prof_cas(nlev_cas))
286	allocate(vq_prof_cas(nlev_cas))
287	allocate(vu_prof_cas(nlev_cas))
288	allocate(vv_prof_cas(nlev_cas))
289	allocate(dt_prof_cas(nlev_cas))
290	allocate(dth_prof_cas(nlev_cas))
291	allocate(dtrad_prof_cas(nlev_cas))
292	allocate(dq_prof_cas(nlev_cas))
293	allocate(du_prof_cas(nlev_cas))
294	allocate(dv_prof_cas(nlev_cas))
295	allocate(uw_prof_cas(nlev_cas))
296	allocate(vw_prof_cas(nlev_cas))
297	allocate(q1_prof_cas(nlev_cas))
298	allocate(q2_prof_cas(nlev_cas))
299
300	print*,'Allocations OK'
301	call read2_cas (nid,nlev_cas,nt_cas, &
302	& ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas, &
303	& ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,ug_cas,vg_cas,du_cas,hu_cas,vu_cas, &
304	& dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas, &
305	& dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tke_cas, &
306	& uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough, &
307	& o3_cas,rugos_cas,clay_cas,sand_cas)
308	print*,'Read2 cas OK'
309	do ii=1,nlev_cas
310	print*,'apres read2_cas, plev_cas=',ii,plev_cas(ii,1)
311	enddo
312
313
314	END SUBROUTINE read2_1D_cas
315
316
317
318	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
319	SUBROUTINE deallocate2_1D_cases
320	!profils environnementaux:
321	deallocate(plev_cas,plevh_cas)
322
323	deallocate(z_cas,zh_cas)
324	deallocate(ap_cas,bp_cas)
325	deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
326	deallocate(th_cas,thl_cas,thv_cas,rv_cas)
327	deallocate(u_cas,v_cas,vitw_cas,omega_cas)
328
329	!forcing
330	deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
331	deallocate(hq_cas,vq_cas,dq_cas)
332	deallocate(hth_cas,vth_cas,dth_cas)
333	deallocate(hr_cas,vr_cas,dr_cas)
334	deallocate(hu_cas,vu_cas,du_cas)
335	deallocate(hv_cas,vv_cas,dv_cas)
336	deallocate(ug_cas)
337	deallocate(vg_cas)
338	deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tke_cas,uw_cas,vw_cas,q1_cas,q2_cas)
339
340	!champs interpoles
341	deallocate(plev_prof_cas)
342	deallocate(t_prof_cas)
343	deallocate(theta_prof_cas)
344	deallocate(thl_prof_cas)
345	deallocate(thv_prof_cas)
346	deallocate(q_prof_cas)
347	deallocate(qv_prof_cas)
348	deallocate(ql_prof_cas)
349	deallocate(qi_prof_cas)
350	deallocate(rh_prof_cas)
351	deallocate(rv_prof_cas)
352	deallocate(u_prof_cas)
353	deallocate(v_prof_cas)
354	deallocate(vitw_prof_cas)
355	deallocate(omega_prof_cas)
356	deallocate(ug_prof_cas)
357	deallocate(vg_prof_cas)
358	deallocate(ht_prof_cas)
359	deallocate(hq_prof_cas)
360	deallocate(hu_prof_cas)
361	deallocate(hv_prof_cas)
362	deallocate(vt_prof_cas)
363	deallocate(vq_prof_cas)
364	deallocate(vu_prof_cas)
365	deallocate(vv_prof_cas)
366	deallocate(dt_prof_cas)
367	deallocate(dtrad_prof_cas)
368	deallocate(dq_prof_cas)
369	deallocate(du_prof_cas)
370	deallocate(dv_prof_cas)
371	deallocate(t_prof_cas)
372	deallocate(u_prof_cas)
373	deallocate(v_prof_cas)
374	deallocate(uw_prof_cas)
375	deallocate(vw_prof_cas)
376	deallocate(q1_prof_cas)
377	deallocate(q2_prof_cas)
378
379	END SUBROUTINE deallocate2_1D_cases
380
381
382	END MODULE mod_1D_cases_read2
383	!=====================================================================
384	subroutine read_cas(nid,nlevel,ntime &
385	& ,zz,pp,temp,qv,rh,theta,rv,u,v,ug,vg,w, &
386	& du,hu,vu,dv,hv,vv,dt,dtrad,ht,vt,dq,hq,vq, &
387	& dth,hth,vth,dr,hr,vr,sens,flat,ts,ustar,uw,vw,q1,q2)
388
389	!program reading forcing of the case study
390	implicit none
391	#include "netcdf.inc"
392
393	integer ntime,nlevel
394
395	real zz(nlevel,ntime)
396	real pp(nlevel,ntime)
397	real temp(nlevel,ntime),qv(nlevel,ntime),rh(nlevel,ntime)
398	real theta(nlevel,ntime),rv(nlevel,ntime)
399	real u(nlevel,ntime)
400	real v(nlevel,ntime)
401	real ug(nlevel,ntime)
402	real vg(nlevel,ntime)
403	real w(nlevel,ntime)
404	real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
405	real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
406	real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
407	real dtrad(nlevel,ntime)
408	real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
409	real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime)
410	real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
411	real flat(ntime),sens(ntime),ts(ntime),ustar(ntime)
412	real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime),resul(nlevel,ntime),resul1(ntime)
413
414
415	integer nid, ierr, ierr1,ierr2,rid,i
416	integer nbvar3d
417	parameter(nbvar3d=39)
418	integer var3didin(nbvar3d)
419	character*5 name_var(1:nbvar3d)
420	data name_var/'zz','pp','temp','qv','rh','theta','rv','u','v','ug','vg','w','advu','hu','vu',&
421	&'advv','hv','vv','advT','hT','vT','advq','hq','vq','advth','hth','vth','advr','hr','vr',&
422	&'radT','uw','vw','q1','q2','sens','flat','ts','ustar'/
423
424	do i=1,nbvar3d
425	print *,'Dans read_cas, on va lire ',nid,i,name_var(i)
426	enddo
427	do i=1,nbvar3d
428	ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
429	print *,'ierr=',i,ierr,name_var(i),var3didin(i)
430	if(ierr/=NF_NOERR) then
431	print *,'Variable manquante dans cas.nc:',name_var(i)
432	endif
433	enddo
434	do i=1,nbvar3d
435	print *,'Dans read_cas, on va lire ',var3didin(i),name_var(i)
436	if(i.LE.35) then
437	#ifdef NC_DOUBLE
438	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
439	#else
440	ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul)
441	#endif
442	print *,'Dans read_cas, on a lu ',ierr,var3didin(i),name_var(i)
443	if(ierr/=NF_NOERR) then
444	print *,'Pb a la lecture de cas.nc: ',name_var(i)
445	stop "getvarup"
446	endif
447	else
448	#ifdef NC_DOUBLE
449	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul1)
450	#else
451	ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul1)
452	#endif
453	print *,'Dans read_cas, on a lu ',ierr,var3didin(i),name_var(i)
454	if(ierr/=NF_NOERR) then
455	print *,'Pb a la lecture de cas.nc: ',name_var(i)
456	stop "getvarup"
457	endif
458	endif
459	select case(i)
460	case(1) ; zz=resul
461	case(2) ; pp=resul
462	case(3) ; temp=resul
463	case(4) ; qv=resul
464	case(5) ; rh=resul
465	case(6) ; theta=resul
466	case(7) ; rv=resul
467	case(8) ; u=resul
468	case(9) ; v=resul
469	case(10) ; ug=resul
470	case(11) ; vg=resul
471	case(12) ; w=resul
472	case(13) ; du=resul
473	case(14) ; hu=resul
474	case(15) ; vu=resul
475	case(16) ; dv=resul
476	case(17) ; hv=resul
477	case(18) ; vv=resul
478	case(19) ; dt=resul
479	case(20) ; ht=resul
480	case(21) ; vt=resul
481	case(22) ; dq=resul
482	case(23) ; hq=resul
483	case(24) ; vq=resul
484	case(25) ; dth=resul
485	case(26) ; hth=resul
486	case(27) ; vth=resul
487	case(28) ; dr=resul
488	case(29) ; hr=resul
489	case(30) ; vr=resul
490	case(31) ; dtrad=resul
491	case(32) ; uw=resul
492	case(33) ; vw=resul
493	case(34) ; q1=resul
494	case(35) ; q2=resul
495	case(36) ; sens=resul1
496	case(37) ; flat=resul1
497	case(38) ; ts=resul1
498	case(39) ; ustar=resul1
499	end select
500	enddo
501
502	return
503	end subroutine read_cas
504	!======================================================================
505	subroutine read2_cas(nid,nlevel,ntime, &
506	& ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
507	& du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq, &
508	& dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tke,uw,vw,q1,q2, &
509	& orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough, &
510	& heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
511
512	!program reading forcing of the case study
513	implicit none
514	#include "netcdf.inc"
515
516	integer ntime,nlevel
517
518	real ap(nlevel+1),bp(nlevel+1)
519	real zz(nlevel,ntime),zzh(nlevel+1)
520	real pp(nlevel,ntime),pph(nlevel+1)
521	real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
522	real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
523	real u(nlevel,ntime),v(nlevel,ntime)
524	real ug(nlevel,ntime),vg(nlevel,ntime)
525	real vitw(nlevel,ntime),omega(nlevel,ntime)
526	real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
527	real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
528	real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
529	real dtrad(nlevel,ntime)
530	real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
531	real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
532	real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
533	real flat(ntime),sens(ntime),ustar(ntime)
534	real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
535	real ts(ntime),ps(ntime),tke(ntime)
536	real orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas
537	real apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
538
539
540	integer nid, ierr,ierr1,ierr2,rid,i
541	integer nbvar3d
542	parameter(nbvar3d=62)
543	integer var3didin(nbvar3d),missing_var(nbvar3d)
544	character*12 name_var(1:nbvar3d)
545	data name_var/'coor_par_a','coor_par_b','height_h','pressure_h',&
546	&'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','tadv','tadvh','tadvv',&
547	&'qadv','qadvh','qadvv','thadv','thadvh','thadvv','thladvh','radv','radvh','radvv','radcool','q1','q2','ustress','vstress','rh',&
548	&'height_f','pressure_f','temp','theta','thv','thl','qv','ql','qi','rv','u','v',&
549	&'sfc_sens_flx','sfc_lat_flx','ts','ps','ustar','tke',&
550	&'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough','o3','rugos','clay','sand'/
551	do i=1,nbvar3d
552	missing_var(i)=0.
553	enddo
554
555	!-----------------------------------------------------------------------
556	do i=1,nbvar3d
557	ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
558	if(ierr/=NF_NOERR) then
559	print *,'Variable manquante dans cas.nc:',i,name_var(i)
560	ierr=NF_NOERR
561	missing_var(i)=1
562	else
563	!-----------------------------------------------------------------------
564	if(i.LE.4) then ! Lecture des coord pression en (nlevelp1,lat,lon)
565	#ifdef NC_DOUBLE
566	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),apbp)
567	#else
568	ierr = NF_GET_VAR_REAL(nid,var3didin(i),apbp)
569	#endif
570	print *,'read2_cas(apbp), on a lu ',i,name_var(i)
571	if(ierr/=NF_NOERR) then
572	print *,'Pb a la lecture de cas.nc: ',name_var(i)
573	stop "getvarup"
574	endif
575	!-----------------------------------------------------------------------
576	else if(i.gt.4.and.i.LE.45) then ! Lecture des variables en (time,nlevel,lat,lon)
577	#ifdef NC_DOUBLE
578	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
579	#else
580	ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul)
581	#endif
582	print *,'read2_cas(resul), on a lu ',i,name_var(i)
583	if(ierr/=NF_NOERR) then
584	print *,'Pb a la lecture de cas.nc: ',name_var(i)
585	stop "getvarup"
586	endif
587	!-----------------------------------------------------------------------
588	else if (i.gt.45.and.i.LE.51) then ! Lecture des variables en (time,lat,lon)
589	#ifdef NC_DOUBLE
590	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul2)
591	#else
592	ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul2)
593	#endif
594	print *,'read2_cas(resul2), on a lu ',i,name_var(i)
595	if(ierr/=NF_NOERR) then
596	print *,'Pb a la lecture de cas.nc: ',name_var(i)
597	stop "getvarup"
598	endif
599	!-----------------------------------------------------------------------
600	else ! Lecture des constantes (lat,lon)
601	#ifdef NC_DOUBLE
602	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul3)
603	#else
604	ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul3)
605	#endif
606	print *,'read2_cas(resul3), on a lu ',i,name_var(i)
607	if(ierr/=NF_NOERR) then
608	print *,'Pb a la lecture de cas.nc: ',name_var(i)
609	stop "getvarup"
610	endif
611	endif
612	endif
613	!-----------------------------------------------------------------------
614	select case(i)
615	case(1) ; ap=apbp ! donnees indexees en nlevel+1
616	case(2) ; bp=apbp
617	case(3) ; zzh=apbp
618	case(4) ; pph=apbp
619	case(5) ; vitw=resul ! donnees indexees en nlevel,time
620	case(6) ; omega=resul
621	case(7) ; ug=resul
622	case(8) ; vg=resul
623	case(9) ; du=resul
624	case(10) ; hu=resul
625	case(11) ; vu=resul
626	case(12) ; dv=resul
627	case(13) ; hv=resul
628	case(14) ; vv=resul
629	case(15) ; dt=resul
630	case(16) ; ht=resul
631	case(17) ; vt=resul
632	case(18) ; dq=resul
633	case(19) ; hq=resul
634	case(20) ; vq=resul
635	case(21) ; dth=resul
636	case(22) ; hth=resul
637	case(23) ; vth=resul
638	case(24) ; hthl=resul
639	case(25) ; dr=resul
640	case(26) ; hr=resul
641	case(27) ; vr=resul
642	case(28) ; dtrad=resul
643	case(29) ; q1=resul
644	case(30) ; q2=resul
645	case(31) ; uw=resul
646	case(32) ; vw=resul
647	case(33) ; rh=resul
648	case(34) ; zz=resul ! donnees en time,nlevel pour profil initial
649	case(35) ; pp=resul
650	case(36) ; temp=resul
651	case(37) ; theta=resul
652	case(38) ; thv=resul
653	case(39) ; thl=resul
654	case(40) ; qv=resul
655	case(41) ; ql=resul
656	case(42) ; qi=resul
657	case(43) ; rv=resul
658	case(44) ; u=resul
659	case(45) ; v=resul
660	case(46) ; sens=resul2 ! donnees indexees en time
661	case(47) ; flat=resul2
662	case(48) ; ts=resul2
663	case(49) ; ps=resul2
664	case(50) ; ustar=resul2
665	case(51) ; tke=resul2
666	case(52) ; orog_cas=resul3 ! constantes
667	case(53) ; albedo_cas=resul3
668	case(54) ; emiss_cas=resul3
669	case(55) ; t_skin_cas=resul3
670	case(56) ; q_skin_cas=resul3
671	case(57) ; mom_rough=resul3
672	case(58) ; heat_rough=resul3
673	case(59) ; o3_cas=resul3
674	case(60) ; rugos_cas=resul3
675	case(61) ; clay_cas=resul3
676	case(62) ; sand_cas=resul3
677	end select
678	resul=0.
679	resul1=0.
680	resul2=0.
681	resul3=0.
682	enddo
683	!-----------------------------------------------------------------------
684
685	return
686	end subroutine read2_cas
687	!======================================================================
688	SUBROUTINE interp_case_time(day,day1,annee_ref &
689	! & ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas &
690	& ,nt_cas,nlev_cas &
691	& ,ts_cas,ps_cas,plev_cas,t_cas,q_cas,u_cas,v_cas &
692	& ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas &
693	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
694	& ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas &
695	& ,uw_cas,vw_cas,q1_cas,q2_cas &
696	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas &
697	& ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas &
698	& ,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas &
699	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas &
700	& ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas &
701	& ,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas &
702	& ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
703
704
705	implicit none
706
707	!---------------------------------------------------------------------------------------
708	! Time interpolation of a 2D field to the timestep corresponding to day
709	!
710	! day: current julian day (e.g. 717538.2)
711	! day1: first day of the simulation
712	! nt_cas: total nb of data in the forcing
713	! pdt_cas: total time interval (in sec) between 2 forcing data
714	!---------------------------------------------------------------------------------------
715
716	#include "compar1d.h"
717	#include "date_cas.h"
718
719	! inputs:
720	integer annee_ref
721	integer nt_cas,nlev_cas
722	real day, day1,day_cas
723	real ts_cas(nt_cas),ps_cas(nt_cas)
724	real plev_cas(nlev_cas,nt_cas)
725	real t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas)
726	real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
727	real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
728	real vitw_cas(nlev_cas,nt_cas)
729	real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
730	real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
731	real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
732	real dtrad_cas(nlev_cas,nt_cas)
733	real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
734	real lat_cas(nt_cas)
735	real sens_cas(nt_cas)
736	real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
737	real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
738
739	! outputs:
740	real plev_prof_cas(nlev_cas)
741	real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
742	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
743	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
744	real vitw_prof_cas(nlev_cas)
745	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
746	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
747	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
748	real dtrad_prof_cas(nlev_cas)
749	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
750	real lat_prof_cas,sens_prof_cas,ts_prof_cas,ustar_prof_cas
751	real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
752	! local:
753	integer it_cas1, it_cas2,k
754	real timeit,time_cas1,time_cas2,frac
755
756
757	print*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
758
759	! On teste si la date du cas AMMA est correcte.
760	! C est pour memoire car en fait les fichiers .def
761	! sont censes etre corrects.
762	! A supprimer a terme (MPL 20150623)
763	! if ((forcing_type.eq.10).and.(1.eq.0)) then
764	! Check that initial day of the simulation consistent with AMMA case:
765	! if (annee_ref.ne.2006) then
766	! print*,'Pour AMMA, annee_ref doit etre 2006'
767	! print*,'Changer annee_ref dans run.def'
768	! stop
769	! endif
770	! if (annee_ref.eq.2006 .and. day1.lt.day_cas) then
771	! print*,'AMMA a debute le 10 juillet 2006',day1,day_cas
772	! print*,'Changer dayref dans run.def'
773	! stop
774	! endif
775	! if (annee_ref.eq.2006 .and. day1.gt.day_cas+1) then
776	! print*,'AMMA a fini le 11 juillet'
777	! print*,'Changer dayref ou nday dans run.def'
778	! stop
779	! endif
780	! endif
781
782	! Determine timestep relative to the 1st day:
783	! timeit=(day-day1)*86400.
784	! if (annee_ref.eq.1992) then
785	! timeit=(day-day_cas)*86400.
786	! else
787	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
788	! endif
789	timeit=(day-day_ju_ini_cas)*86400
790	print *,'day=',day
791	print *,'day_ju_ini_cas=',day_ju_ini_cas
792	print *,'pdt_cas=',pdt_cas
793	print *,'timeit=',timeit
794	print *,'nt_cas=',nt_cas
795
796	! Determine the closest observation times:
797	! it_cas1=INT(timeit/pdt_cas)+1
798	! it_cas2=it_cas1 + 1
799	! time_cas1=(it_cas1-1)*pdt_cas
800	! time_cas2=(it_cas2-1)*pdt_cas
801
802	it_cas1=INT(timeit/pdt_cas)+1
803	IF (it_cas1 .EQ. nt_cas) THEN
804	it_cas2=it_cas1
805	ELSE
806	it_cas2=it_cas1 + 1
807	ENDIF
808	time_cas1=(it_cas1-1)*pdt_cas
809	time_cas2=(it_cas2-1)*pdt_cas
810	print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
811
812	if (it_cas1 .gt. nt_cas) then
813	write(,) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: ' &
814	& ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
815	stop
816	endif
817
818	! time interpolation:
819	IF (it_cas1 .EQ. it_cas2) THEN
820	frac=0.
821	ELSE
822	frac=(time_cas2-timeit)/(time_cas2-time_cas1)
823	frac=max(frac,0.0)
824	ENDIF
825
826	lat_prof_cas = lat_cas(it_cas2) &
827	& -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
828	sens_prof_cas = sens_cas(it_cas2) &
829	& -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
830	ts_prof_cas = ts_cas(it_cas2) &
831	& -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
832	ustar_prof_cas = ustar_cas(it_cas2) &
833	& -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
834
835	do k=1,nlev_cas
836	plev_prof_cas(k) = plev_cas(k,it_cas2) &
837	& -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
838	t_prof_cas(k) = t_cas(k,it_cas2) &
839	& -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
840	q_prof_cas(k) = q_cas(k,it_cas2) &
841	& -frac*(q_cas(k,it_cas2)-q_cas(k,it_cas1))
842	u_prof_cas(k) = u_cas(k,it_cas2) &
843	& -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
844	v_prof_cas(k) = v_cas(k,it_cas2) &
845	& -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
846	ug_prof_cas(k) = ug_cas(k,it_cas2) &
847	& -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
848	vg_prof_cas(k) = vg_cas(k,it_cas2) &
849	& -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
850	vitw_prof_cas(k) = vitw_cas(k,it_cas2) &
851	& -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
852	du_prof_cas(k) = du_cas(k,it_cas2) &
853	& -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
854	hu_prof_cas(k) = hu_cas(k,it_cas2) &
855	& -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
856	vu_prof_cas(k) = vu_cas(k,it_cas2) &
857	& -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
858	dv_prof_cas(k) = dv_cas(k,it_cas2) &
859	& -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
860	hv_prof_cas(k) = hv_cas(k,it_cas2) &
861	& -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
862	vv_prof_cas(k) = vv_cas(k,it_cas2) &
863	& -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
864	dt_prof_cas(k) = dt_cas(k,it_cas2) &
865	& -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
866	ht_prof_cas(k) = ht_cas(k,it_cas2) &
867	& -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
868	vt_prof_cas(k) = vt_cas(k,it_cas2) &
869	& -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
870	dtrad_prof_cas(k) = dtrad_cas(k,it_cas2) &
871	& -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
872	dq_prof_cas(k) = dq_cas(k,it_cas2) &
873	& -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
874	hq_prof_cas(k) = hq_cas(k,it_cas2) &
875	& -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
876	vq_prof_cas(k) = vq_cas(k,it_cas2) &
877	& -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
878	uw_prof_cas(k) = uw_cas(k,it_cas2) &
879	& -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
880	vw_prof_cas(k) = vw_cas(k,it_cas2) &
881	& -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
882	q1_prof_cas(k) = q1_cas(k,it_cas2) &
883	& -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
884	q2_prof_cas(k) = q2_cas(k,it_cas2) &
885	& -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
886	enddo
887
888	return
889	END SUBROUTINE interp_case_time
890
891	!**********************************************************************************************
892	SUBROUTINE interp2_case_time(day,day1,annee_ref &
893	! & ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas &
894	& ,nt_cas,nlev_cas &
895	& ,ts_cas,ps_cas,plev_cas,t_cas,theta_cas,thv_cas,thl_cas &
896	& ,qv_cas,ql_cas,qi_cas,u_cas,v_cas &
897	& ,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas &
898	& ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas &
899	& ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas &
900	& ,lat_cas,sens_cas,ustar_cas &
901	& ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas &
902	!
903	& ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas &
904	& ,thv_prof_cas,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas &
905	& ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas &
906	& ,vitw_prof_cas,omega_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas &
907	& ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas &
908	& ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas &
909	& ,hq_prof_cas,vq_prof_cas,dth_prof_cas,hth_prof_cas,vth_prof_cas &
910	& ,lat_prof_cas,sens_prof_cas &
911	& ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
912
913
914	implicit none
915
916	!---------------------------------------------------------------------------------------
917	! Time interpolation of a 2D field to the timestep corresponding to day
918	!
919	! day: current julian day (e.g. 717538.2)
920	! day1: first day of the simulation
921	! nt_cas: total nb of data in the forcing
922	! pdt_cas: total time interval (in sec) between 2 forcing data
923	!---------------------------------------------------------------------------------------
924
925	#include "compar1d.h"
926	#include "date_cas.h"
927
928	! inputs:
929	integer annee_ref
930	integer nt_cas,nlev_cas
931	real day, day1,day_cas
932	real ts_cas(nt_cas),ps_cas(nt_cas)
933	real plev_cas(nlev_cas,nt_cas)
934	real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas)
935	real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
936	real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
937	real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
938	real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas)
939	real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
940	real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
941	real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
942	real dth_cas(nlev_cas,nt_cas),hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas)
943	real dtrad_cas(nlev_cas,nt_cas)
944	real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
945	real lat_cas(nt_cas),sens_cas(nt_cas),tke_cas(nt_cas)
946	real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
947	real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
948
949	! outputs:
950	real plev_prof_cas(nlev_cas)
951	real t_prof_cas(nlev_cas),theta_prof_cas(nlev_cas),thl_prof_cas(nlev_cas),thv_prof_cas(nlev_cas)
952	real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
953	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
954	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
955	real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
956	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
957	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
958	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
959	real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
960	real dtrad_prof_cas(nlev_cas)
961	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
962	real lat_prof_cas,sens_prof_cas,tke_prof_cas,ts_prof_cas,ustar_prof_cas
963	real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
964	! local:
965	integer it_cas1, it_cas2,k
966	real timeit,time_cas1,time_cas2,frac
967
968
969	print*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
970	! do k=1,nlev_cas
971	! print*,'debut de interp2_case_time, plev_cas=',k,plev_cas(k,1)
972	! enddo
973
974	! On teste si la date du cas AMMA est correcte.
975	! C est pour memoire car en fait les fichiers .def
976	! sont censes etre corrects.
977	! A supprimer a terme (MPL 20150623)
978	! if ((forcing_type.eq.10).and.(1.eq.0)) then
979	! Check that initial day of the simulation consistent with AMMA case:
980	! if (annee_ref.ne.2006) then
981	! print*,'Pour AMMA, annee_ref doit etre 2006'
982	! print*,'Changer annee_ref dans run.def'
983	! stop
984	! endif
985	! if (annee_ref.eq.2006 .and. day1.lt.day_cas) then
986	! print*,'AMMA a debute le 10 juillet 2006',day1,day_cas
987	! print*,'Changer dayref dans run.def'
988	! stop
989	! endif
990	! if (annee_ref.eq.2006 .and. day1.gt.day_cas+1) then
991	! print*,'AMMA a fini le 11 juillet'
992	! print*,'Changer dayref ou nday dans run.def'
993	! stop
994	! endif
995	! endif
996
997	! Determine timestep relative to the 1st day:
998	! timeit=(day-day1)*86400.
999	! if (annee_ref.eq.1992) then
1000	! timeit=(day-day_cas)*86400.
1001	! else
1002	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
1003	! endif
1004	timeit=(day-day_ju_ini_cas)*86400
1005	print *,'day=',day
1006	print *,'day_ju_ini_cas=',day_ju_ini_cas
1007	print *,'pdt_cas=',pdt_cas
1008	print *,'timeit=',timeit
1009	print *,'nt_cas=',nt_cas
1010
1011	! Determine the closest observation times:
1012	! it_cas1=INT(timeit/pdt_cas)+1
1013	! it_cas2=it_cas1 + 1
1014	! time_cas1=(it_cas1-1)*pdt_cas
1015	! time_cas2=(it_cas2-1)*pdt_cas
1016
1017	it_cas1=INT(timeit/pdt_cas)+1
1018	IF (it_cas1 .EQ. nt_cas) THEN
1019	it_cas2=it_cas1
1020	ELSE
1021	it_cas2=it_cas1 + 1
1022	ENDIF
1023	time_cas1=(it_cas1-1)*pdt_cas
1024	time_cas2=(it_cas2-1)*pdt_cas
1025	print *,'timeit,pdt_cas,nt_cas=',timeit,pdt_cas,nt_cas
1026	print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
1027
1028	if (it_cas1 .gt. nt_cas) then
1029	write(,) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: ' &
1030	& ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
1031	stop
1032	endif
1033
1034	! time interpolation:
1035	IF (it_cas1 .EQ. it_cas2) THEN
1036	frac=0.
1037	ELSE
1038	frac=(time_cas2-timeit)/(time_cas2-time_cas1)
1039	frac=max(frac,0.0)
1040	ENDIF
1041
1042	lat_prof_cas = lat_cas(it_cas2) &
1043	& -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
1044	sens_prof_cas = sens_cas(it_cas2) &
1045	& -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
1046	tke_prof_cas = tke_cas(it_cas2) &
1047	& -frac*(tke_cas(it_cas2)-tke_cas(it_cas1))
1048	ts_prof_cas = ts_cas(it_cas2) &
1049	& -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
1050	ustar_prof_cas = ustar_cas(it_cas2) &
1051	& -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
1052
1053	do k=1,nlev_cas
1054	plev_prof_cas(k) = plev_cas(k,it_cas2) &
1055	& -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
1056	t_prof_cas(k) = t_cas(k,it_cas2) &
1057	& -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
1058	print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
1059	theta_prof_cas(k) = theta_cas(k,it_cas2) &
1060	& -frac*(theta_cas(k,it_cas2)-theta_cas(k,it_cas1))
1061	thv_prof_cas(k) = thv_cas(k,it_cas2) &
1062	& -frac*(thv_cas(k,it_cas2)-thv_cas(k,it_cas1))
1063	thl_prof_cas(k) = thl_cas(k,it_cas2) &
1064	& -frac*(thl_cas(k,it_cas2)-thl_cas(k,it_cas1))
1065	qv_prof_cas(k) = qv_cas(k,it_cas2) &
1066	& -frac*(qv_cas(k,it_cas2)-qv_cas(k,it_cas1))
1067	ql_prof_cas(k) = ql_cas(k,it_cas2) &
1068	& -frac*(ql_cas(k,it_cas2)-ql_cas(k,it_cas1))
1069	qi_prof_cas(k) = qi_cas(k,it_cas2) &
1070	& -frac*(qi_cas(k,it_cas2)-qi_cas(k,it_cas1))
1071	u_prof_cas(k) = u_cas(k,it_cas2) &
1072	& -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
1073	v_prof_cas(k) = v_cas(k,it_cas2) &
1074	& -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
1075	ug_prof_cas(k) = ug_cas(k,it_cas2) &
1076	& -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
1077	vg_prof_cas(k) = vg_cas(k,it_cas2) &
1078	& -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
1079	vitw_prof_cas(k) = vitw_cas(k,it_cas2) &
1080	& -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
1081	omega_prof_cas(k) = omega_cas(k,it_cas2) &
1082	& -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
1083	du_prof_cas(k) = du_cas(k,it_cas2) &
1084	& -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
1085	hu_prof_cas(k) = hu_cas(k,it_cas2) &
1086	& -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
1087	vu_prof_cas(k) = vu_cas(k,it_cas2) &
1088	& -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
1089	dv_prof_cas(k) = dv_cas(k,it_cas2) &
1090	& -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
1091	hv_prof_cas(k) = hv_cas(k,it_cas2) &
1092	& -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
1093	vv_prof_cas(k) = vv_cas(k,it_cas2) &
1094	& -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
1095	dt_prof_cas(k) = dt_cas(k,it_cas2) &
1096	& -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
1097	ht_prof_cas(k) = ht_cas(k,it_cas2) &
1098	& -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
1099	vt_prof_cas(k) = vt_cas(k,it_cas2) &
1100	& -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
1101	dth_prof_cas(k) = dth_cas(k,it_cas2) &
1102	& -frac*(dth_cas(k,it_cas2)-dth_cas(k,it_cas1))
1103	hth_prof_cas(k) = hth_cas(k,it_cas2) &
1104	& -frac*(hth_cas(k,it_cas2)-hth_cas(k,it_cas1))
1105	vth_prof_cas(k) = vth_cas(k,it_cas2) &
1106	& -frac*(vth_cas(k,it_cas2)-vth_cas(k,it_cas1))
1107	dtrad_prof_cas(k) = dtrad_cas(k,it_cas2) &
1108	& -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
1109	dq_prof_cas(k) = dq_cas(k,it_cas2) &
1110	& -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
1111	hq_prof_cas(k) = hq_cas(k,it_cas2) &
1112	& -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
1113	vq_prof_cas(k) = vq_cas(k,it_cas2) &
1114	& -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
1115	uw_prof_cas(k) = uw_cas(k,it_cas2) &
1116	& -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
1117	vw_prof_cas(k) = vw_cas(k,it_cas2) &
1118	& -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
1119	q1_prof_cas(k) = q1_cas(k,it_cas2) &
1120	& -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
1121	q2_prof_cas(k) = q2_cas(k,it_cas2) &
1122	& -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
1123	enddo
1124
1125	return
1126	END SUBROUTINE interp2_case_time
1127
1128	!**********************************************************************************************
1129

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

Download in other formats:

Original Format