Context Navigation

mod_1D_cases_read2.F90 @ 5403

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

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

Download in other formats:

Original Format