Context Navigation

source: LMDZ6/trunk/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h @ 4824

Last change on this file since 4824 was 4593, checked in by yann meurdesoif, 18 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

File size: 108.1 KB

Line
1	!======================================================================
2	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &
3	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
4	& ,ht_toga,vt_toga,hq_toga,vq_toga)
5	implicit none
6
7	!-------------------------------------------------------------------------
8	! Read TOGA-COARE forcing data
9	!-------------------------------------------------------------------------
10
11	integer nlev_toga,nt_toga
12	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
13	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
14	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
15	real w_toga(nlev_toga,nt_toga)
16	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
17	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
18	character*80 fich_toga
19
20	integer k,ip
21	real bid
22
23	integer iy,im,id,ih
24
25	real plev_min
26
27	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
28
29	open(21,file=trim(fich_toga),form='formatted')
30	read(21,'(a)')
31	do ip = 1, nt_toga
32	read(21,'(a)')
33	read(21,'(a)')
34	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
35	read(21,'(a)')
36	read(21,'(a)')
37
38	do k = 1, nlev_toga
39	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &
40	& ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &
41	& ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
42
43	! conversion in SI units:
44	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
45	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
46	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
47	! no water vapour tendency above 55 hPa
48	if (plev_toga(k,ip) .lt. plev_min) then
49	q_toga(k,ip) = 0.
50	hq_toga(k,ip) = 0.
51	vq_toga(k,ip) =0.
52	endif
53	enddo
54
55	ts_toga(ip)=ts_toga(ip)+273.15 ! K
56	enddo
57	close(21)
58
59	223 format(4i3,6f8.2)
60	230 format(6f9.3,4e11.3)
61
62	return
63	end
64
65	!-------------------------------------------------------------------------
66	SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
67	implicit none
68
69	!-------------------------------------------------------------------------
70	! Read I.SANDU case forcing data
71	!-------------------------------------------------------------------------
72
73	integer nlev_sandu,nt_sandu
74	real ts_sandu(nt_sandu)
75	character*80 fich_sandu
76
77	integer ip
78	integer iy,im,id,ih
79
80	real plev_min
81
82	print*,'nlev_sandu',nlev_sandu
83	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
84
85	open(21,file=trim(fich_sandu),form='formatted')
86	read(21,'(a)')
87	do ip = 1, nt_sandu
88	read(21,'(a)')
89	read(21,'(a)')
90	read(21,223) iy, im, id, ih, ts_sandu(ip)
91	print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
92	enddo
93	close(21)
94
95	223 format(4i3,f8.2)
96
97	return
98	end
99
100	!=====================================================================
101	!-------------------------------------------------------------------------
102	SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, &
103	& ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
104	implicit none
105
106	!-------------------------------------------------------------------------
107	! Read Astex case forcing data
108	!-------------------------------------------------------------------------
109
110	integer nlev_astex,nt_astex
111	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
112	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
113	character*80 fich_astex
114
115	integer ip
116	integer iy,im,id,ih
117
118	real plev_min
119
120	print*,'nlev_astex',nlev_astex
121	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
122
123	open(21,file=trim(fich_astex),form='formatted')
124	read(21,'(a)')
125	read(21,'(a)')
126	do ip = 1, nt_astex
127	read(21,'(a)')
128	read(21,'(a)')
129	read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), &
130	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
131	ts_astex(ip)=ts_astex(ip)+273.15
132	print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), &
133	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
134	enddo
135	close(21)
136
137	223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
138
139	return
140	end
141	!=====================================================================
142	subroutine read_twpice(fich_twpice,nlevel,ntime &
143	& ,T_srf,plev,T,q,u,v,omega &
144	& ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
145
146	!program reading forcings of the TWP-ICE experiment
147
148	use netcdf, only: nf90_get_var
149
150	implicit none
151
152	INCLUDE "netcdf.inc"
153
154	integer ntime,nlevel
155	integer l,k
156	character*80 :: fich_twpice
157	real*8 time(ntime)
158	real*8 lat, lon, alt, phis
159	real*8 lev(nlevel)
160	real*8 plev(nlevel,ntime)
161
162	real*8 T(nlevel,ntime)
163	real*8 q(nlevel,ntime),u(nlevel,ntime)
164	real*8 v(nlevel,ntime)
165	real*8 omega(nlevel,ntime), div(nlevel,ntime)
166	real*8 T_adv_h(nlevel,ntime)
167	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
168	real*8 q_adv_v(nlevel,ntime)
169	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
170	real*8 s_adv_v(nlevel,ntime)
171	real*8 p_srf_aver(ntime), p_srf_center(ntime)
172	real*8 T_srf(ntime)
173
174	integer nid, ierr
175	integer nbvar3d
176	parameter(nbvar3d=20)
177	integer var3didin(nbvar3d)
178
179	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
180	if (ierr.NE.NF_NOERR) then
181	write(,) 'ERROR: Pb opening forcings cdf file '
182	write(,) NF_STRERROR(ierr)
183	stop ""
184	endif
185
186	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
187	if(ierr/=NF_NOERR) then
188	write(,) NF_STRERROR(ierr)
189	stop 'lat'
190	endif
191
192	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
193	if(ierr/=NF_NOERR) then
194	write(,) NF_STRERROR(ierr)
195	stop 'lon'
196	endif
197
198	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
199	if(ierr/=NF_NOERR) then
200	write(,) NF_STRERROR(ierr)
201	stop 'alt'
202	endif
203
204	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
205	if(ierr/=NF_NOERR) then
206	write(,) NF_STRERROR(ierr)
207	stop 'phis'
208	endif
209
210	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
211	if(ierr/=NF_NOERR) then
212	write(,) NF_STRERROR(ierr)
213	stop 'T'
214	endif
215
216	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
217	if(ierr/=NF_NOERR) then
218	write(,) NF_STRERROR(ierr)
219	stop 'q'
220	endif
221
222	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
223	if(ierr/=NF_NOERR) then
224	write(,) NF_STRERROR(ierr)
225	stop 'u'
226	endif
227
228	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
229	if(ierr/=NF_NOERR) then
230	write(,) NF_STRERROR(ierr)
231	stop 'v'
232	endif
233
234	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
235	if(ierr/=NF_NOERR) then
236	write(,) NF_STRERROR(ierr)
237	stop 'omega'
238	endif
239
240	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
241	if(ierr/=NF_NOERR) then
242	write(,) NF_STRERROR(ierr)
243	stop 'div'
244	endif
245
246	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
247	if(ierr/=NF_NOERR) then
248	write(,) NF_STRERROR(ierr)
249	stop 'T_adv_h'
250	endif
251
252	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
253	if(ierr/=NF_NOERR) then
254	write(,) NF_STRERROR(ierr)
255	stop 'T_adv_v'
256	endif
257
258	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
259	if(ierr/=NF_NOERR) then
260	write(,) NF_STRERROR(ierr)
261	stop 'q_adv_h'
262	endif
263
264	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
265	if(ierr/=NF_NOERR) then
266	write(,) NF_STRERROR(ierr)
267	stop 'q_adv_v'
268	endif
269
270	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
271	if(ierr/=NF_NOERR) then
272	write(,) NF_STRERROR(ierr)
273	stop 's'
274	endif
275
276	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
277	if(ierr/=NF_NOERR) then
278	write(,) NF_STRERROR(ierr)
279	stop 's_adv_h'
280	endif
281
282	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
283	if(ierr/=NF_NOERR) then
284	write(,) NF_STRERROR(ierr)
285	stop 's_adv_v'
286	endif
287
288	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
289	if(ierr/=NF_NOERR) then
290	write(,) NF_STRERROR(ierr)
291	stop 'p_srf_aver'
292	endif
293
294	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
295	if(ierr/=NF_NOERR) then
296	write(,) NF_STRERROR(ierr)
297	stop 'p_srf_center'
298	endif
299
300	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
301	if(ierr/=NF_NOERR) then
302	write(,) NF_STRERROR(ierr)
303	stop 'T_srf'
304	endif
305
306	!dimensions lecture
307	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
308
309	!pressure
310	do l=1,ntime
311	do k=1,nlevel
312	plev(k,l)=lev(k)
313	enddo
314	enddo
315
316	ierr = NF90_GET_VAR(nid,var3didin(1),lat)
317	if(ierr/=NF_NOERR) then
318	write(,) NF_STRERROR(ierr)
319	stop "getvarup"
320	endif
321	! write(,)'lecture lat ok',lat
322
323	ierr = NF90_GET_VAR(nid,var3didin(2),lon)
324	if(ierr/=NF_NOERR) then
325	write(,) NF_STRERROR(ierr)
326	stop "getvarup"
327	endif
328	! write(,)'lecture lon ok',lon
329
330	ierr = NF90_GET_VAR(nid,var3didin(3),alt)
331	if(ierr/=NF_NOERR) then
332	write(,) NF_STRERROR(ierr)
333	stop "getvarup"
334	endif
335	! write(,)'lecture alt ok',alt
336
337	ierr = NF90_GET_VAR(nid,var3didin(4),phis)
338	if(ierr/=NF_NOERR) then
339	write(,) NF_STRERROR(ierr)
340	stop "getvarup"
341	endif
342	! write(,)'lecture phis ok',phis
343
344	ierr = NF90_GET_VAR(nid,var3didin(5),T)
345	if(ierr/=NF_NOERR) then
346	write(,) NF_STRERROR(ierr)
347	stop "getvarup"
348	endif
349	! write(,)'lecture T ok'
350
351	ierr = NF90_GET_VAR(nid,var3didin(6),q)
352	if(ierr/=NF_NOERR) then
353	write(,) NF_STRERROR(ierr)
354	stop "getvarup"
355	endif
356	! write(,)'lecture q ok'
357	!q in kg/kg
358	do l=1,ntime
359	do k=1,nlevel
360	q(k,l)=q(k,l)/1000.
361	enddo
362	enddo
363	ierr = NF90_GET_VAR(nid,var3didin(7),u)
364	if(ierr/=NF_NOERR) then
365	write(,) NF_STRERROR(ierr)
366	stop "getvarup"
367	endif
368	! write(,)'lecture u ok'
369
370	ierr = NF90_GET_VAR(nid,var3didin(8),v)
371	if(ierr/=NF_NOERR) then
372	write(,) NF_STRERROR(ierr)
373	stop "getvarup"
374	endif
375	! write(,)'lecture v ok'
376
377	ierr = NF90_GET_VAR(nid,var3didin(9),omega)
378	if(ierr/=NF_NOERR) then
379	write(,) NF_STRERROR(ierr)
380	stop "getvarup"
381	endif
382	! write(,)'lecture omega ok'
383	!omega in mb/hour
384	do l=1,ntime
385	do k=1,nlevel
386	omega(k,l)=omega(k,l)*100./3600.
387	enddo
388	enddo
389
390	ierr = NF90_GET_VAR(nid,var3didin(10),div)
391	if(ierr/=NF_NOERR) then
392	write(,) NF_STRERROR(ierr)
393	stop "getvarup"
394	endif
395	! write(,)'lecture div ok'
396
397	ierr = NF90_GET_VAR(nid,var3didin(11),T_adv_h)
398	if(ierr/=NF_NOERR) then
399	write(,) NF_STRERROR(ierr)
400	stop "getvarup"
401	endif
402	! write(,)'lecture T_adv_h ok'
403	!T adv in K/s
404	do l=1,ntime
405	do k=1,nlevel
406	T_adv_h(k,l)=T_adv_h(k,l)/3600.
407	enddo
408	enddo
409
410
411	ierr = NF90_GET_VAR(nid,var3didin(12),T_adv_v)
412	if(ierr/=NF_NOERR) then
413	write(,) NF_STRERROR(ierr)
414	stop "getvarup"
415	endif
416	! write(,)'lecture T_adv_v ok'
417	!T adv in K/s
418	do l=1,ntime
419	do k=1,nlevel
420	T_adv_v(k,l)=T_adv_v(k,l)/3600.
421	enddo
422	enddo
423
424	ierr = NF90_GET_VAR(nid,var3didin(13),q_adv_h)
425	if(ierr/=NF_NOERR) then
426	write(,) NF_STRERROR(ierr)
427	stop "getvarup"
428	endif
429	! write(,)'lecture q_adv_h ok'
430	!q adv in kg/kg/s
431	do l=1,ntime
432	do k=1,nlevel
433	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
434	enddo
435	enddo
436
437
438	ierr = NF90_GET_VAR(nid,var3didin(14),q_adv_v)
439	if(ierr/=NF_NOERR) then
440	write(,) NF_STRERROR(ierr)
441	stop "getvarup"
442	endif
443	! write(,)'lecture q_adv_v ok'
444	!q adv in kg/kg/s
445	do l=1,ntime
446	do k=1,nlevel
447	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
448	enddo
449	enddo
450
451
452	ierr = NF90_GET_VAR(nid,var3didin(15),s)
453	if(ierr/=NF_NOERR) then
454	write(,) NF_STRERROR(ierr)
455	stop "getvarup"
456	endif
457
458	ierr = NF90_GET_VAR(nid,var3didin(16),s_adv_h)
459	if(ierr/=NF_NOERR) then
460	write(,) NF_STRERROR(ierr)
461	stop "getvarup"
462	endif
463
464	ierr = NF90_GET_VAR(nid,var3didin(17),s_adv_v)
465	if(ierr/=NF_NOERR) then
466	write(,) NF_STRERROR(ierr)
467	stop "getvarup"
468	endif
469
470	ierr = NF90_GET_VAR(nid,var3didin(18),p_srf_aver)
471	if(ierr/=NF_NOERR) then
472	write(,) NF_STRERROR(ierr)
473	stop "getvarup"
474	endif
475
476	ierr = NF90_GET_VAR(nid,var3didin(19),p_srf_center)
477	if(ierr/=NF_NOERR) then
478	write(,) NF_STRERROR(ierr)
479	stop "getvarup"
480	endif
481
482	ierr = NF90_GET_VAR(nid,var3didin(20),T_srf)
483	if(ierr/=NF_NOERR) then
484	write(,) NF_STRERROR(ierr)
485	stop "getvarup"
486	endif
487	! write(,)'lecture T_srf ok', T_srf
488
489	return
490	end subroutine read_twpice
491	!=====================================================================
492	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
493
494	use netcdf, only: nf90_get_var
495
496	implicit none
497	INCLUDE "netcdf.inc"
498	integer nid,ttm,llm
499	real*8 time(ttm)
500	real*8 lev(llm)
501	integer ierr
502
503	integer timevar,levvar
504	integer timelen,levlen
505	integer timedimin,levdimin
506
507	! Control & lecture on dimensions
508	! ===============================
509	ierr=NF_INQ_DIMID(nid,"time",timedimin)
510	ierr=NF_INQ_VARID(nid,"time",timevar)
511	if (ierr.NE.NF_NOERR) then
512	write(,) 'ERROR: Field <time> is missing'
513	stop ""
514	endif
515	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
516
517	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
518	ierr=NF_INQ_VARID(nid,"lev",levvar)
519	if (ierr.NE.NF_NOERR) then
520	write(,) 'ERROR: Field <lev> is lacking'
521	stop ""
522	endif
523	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
524
525	if((timelen/=ttm).or.(levlen/=llm)) then
526	write(,) 'ERROR: Not the good lenght for axis'
527	write(,) 'longitude: ',timelen,ttm+1
528	write(,) 'latitude: ',levlen,llm
529	stop ""
530	endif
531
532	ierr = NF90_GET_VAR(nid,timevar,time)
533	ierr = NF90_GET_VAR(nid,levvar,lev)
534
535	return
536	end
537	!=====================================================================
538
539	SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof &
540	& ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof &
541	& ,omega_prof,o3mmr_prof &
542	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
543	& ,omega_mod,o3mmr_mod,mxcalc)
544
545	implicit none
546
547	INCLUDE "dimensions.h"
548
549	!-------------------------------------------------------------------------
550	! Vertical interpolation of SANDUREF forcing data onto model levels
551	!-------------------------------------------------------------------------
552
553	integer nlevmax
554	parameter (nlevmax=41)
555	integer nlev_sandu,mxcalc
556	! real play(llm), plev_prof(nlevmax)
557	! real t_prof(nlevmax),q_prof(nlevmax)
558	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
559	! real ht_prof(nlevmax),vt_prof(nlevmax)
560	! real hq_prof(nlevmax),vq_prof(nlevmax)
561
562	real play(llm), plev_prof(nlev_sandu)
563	real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
564	real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
565	real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
566
567	real t_mod(llm),thl_mod(llm),q_mod(llm)
568	real u_mod(llm),v_mod(llm), w_mod(llm)
569	real omega_mod(llm),o3mmr_mod(llm)
570
571	integer l,k,k1,k2
572	real frac,frac1,frac2,fact
573
574	do l = 1, llm
575
576	if (play(l).ge.plev_prof(nlev_sandu)) then
577
578	mxcalc=l
579	k1=0
580	k2=0
581
582	if (play(l).le.plev_prof(1)) then
583
584	do k = 1, nlev_sandu-1
585	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
586	k1=k
587	k2=k+1
588	endif
589	enddo
590
591	if (k1.eq.0 .or. k2.eq.0) then
592	write(,) 'PB! k1, k2 = ',k1,k2
593	write(,) 'l,play(l) = ',l,play(l)/100
594	do k = 1, nlev_sandu-1
595	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
596	enddo
597	endif
598
599	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
600	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
601	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
602	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
603	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
604	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
605	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
606	omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
607	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
608
609	else !play>plev_prof(1)
610
611	k1=1
612	k2=2
613	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
614	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
615	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
616	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
617	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
618	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
619	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
620	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
621	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
622	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
623
624	endif ! play.le.plev_prof(1)
625
626	else ! above max altitude of forcing file
627
628	!jyg
629	fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
630	fact = max(fact,0.) !jyg
631	fact = exp(-fact) !jyg
632	t_mod(l)= t_prof(nlev_sandu) !jyg
633	thl_mod(l)= thl_prof(nlev_sandu) !jyg
634	q_mod(l)= q_prof(nlev_sandu)*fact !jyg
635	u_mod(l)= u_prof(nlev_sandu)*fact !jyg
636	v_mod(l)= v_prof(nlev_sandu)*fact !jyg
637	w_mod(l)= w_prof(nlev_sandu)*fact !jyg
638	omega_mod(l)= omega_prof(nlev_sandu)*fact !jyg
639	o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact !jyg
640
641	endif ! play
642
643	enddo ! l
644
645	do l = 1,llm
646	! print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
647	! $ l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
648	enddo
649
650	return
651	end
652	!=====================================================================
653	SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof &
654	& ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof &
655	& ,w_prof,tke_prof,o3mmr_prof &
656	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
657	& ,tke_mod,o3mmr_mod,mxcalc)
658
659	implicit none
660
661	INCLUDE "dimensions.h"
662
663	!-------------------------------------------------------------------------
664	! Vertical interpolation of Astex forcing data onto model levels
665	!-------------------------------------------------------------------------
666
667	integer nlevmax
668	parameter (nlevmax=41)
669	integer nlev_astex,mxcalc
670	! real play(llm), plev_prof(nlevmax)
671	! real t_prof(nlevmax),qv_prof(nlevmax)
672	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
673	! real ht_prof(nlevmax),vt_prof(nlevmax)
674	! real hq_prof(nlevmax),vq_prof(nlevmax)
675
676	real play(llm), plev_prof(nlev_astex)
677	real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
678	real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
679	real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
680	real qt_prof(nlev_astex),tke_prof(nlev_astex)
681
682	real t_mod(llm),thl_mod(llm),qv_mod(llm)
683	real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
684	real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
685
686	integer l,k,k1,k2
687	real frac,frac1,frac2,fact
688
689	do l = 1, llm
690
691	if (play(l).ge.plev_prof(nlev_astex)) then
692
693	mxcalc=l
694	k1=0
695	k2=0
696
697	if (play(l).le.plev_prof(1)) then
698
699	do k = 1, nlev_astex-1
700	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
701	k1=k
702	k2=k+1
703	endif
704	enddo
705
706	if (k1.eq.0 .or. k2.eq.0) then
707	write(,) 'PB! k1, k2 = ',k1,k2
708	write(,) 'l,play(l) = ',l,play(l)/100
709	do k = 1, nlev_astex-1
710	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
711	enddo
712	endif
713
714	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
715	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
716	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
717	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
718	ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
719	qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
720	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
721	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
722	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
723	tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
724	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
725
726	else !play>plev_prof(1)
727
728	k1=1
729	k2=2
730	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
731	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
732	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
733	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
734	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
735	ql_mod(l)= frac1ql_prof(k1) - frac2ql_prof(k2)
736	qt_mod(l)= frac1qt_prof(k1) - frac2qt_prof(k2)
737	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
738	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
739	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
740	tke_mod(l)= frac1tke_prof(k1) - frac2tke_prof(k2)
741	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
742
743	endif ! play.le.plev_prof(1)
744
745	else ! above max altitude of forcing file
746
747	!jyg
748	fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
749	fact = max(fact,0.) !jyg
750	fact = exp(-fact) !jyg
751	t_mod(l)= t_prof(nlev_astex) !jyg
752	thl_mod(l)= thl_prof(nlev_astex) !jyg
753	qv_mod(l)= qv_prof(nlev_astex)*fact !jyg
754	ql_mod(l)= ql_prof(nlev_astex)*fact !jyg
755	qt_mod(l)= qt_prof(nlev_astex)*fact !jyg
756	u_mod(l)= u_prof(nlev_astex)*fact !jyg
757	v_mod(l)= v_prof(nlev_astex)*fact !jyg
758	w_mod(l)= w_prof(nlev_astex)*fact !jyg
759	tke_mod(l)= tke_prof(nlev_astex)*fact !jyg
760	o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact !jyg
761
762	endif ! play
763
764	enddo ! l
765
766	do l = 1,llm
767	! print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
768	! $ l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
769	enddo
770
771	return
772	end
773
774	!======================================================================
775	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play &
776	& ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &
777	& ,dth_dyn,dqh_dyn)
778	implicit none
779
780	!-------------------------------------------------------------------------
781	! Read RICO forcing data
782	!-------------------------------------------------------------------------
783	INCLUDE "dimensions.h"
784
785
786	integer nlev_rico
787	real ts_rico,ps_rico
788	real t_rico(llm),q_rico(llm)
789	real u_rico(llm),v_rico(llm)
790	real w_rico(llm)
791	real dth_dyn(llm)
792	real dqh_dyn(llm)
793
794
795	real play(llm),zlay(llm)
796
797
798	real prico(nlev_rico),zrico(nlev_rico)
799
800	character*80 fich_rico
801
802	integer k,l
803
804
805	print*,fich_rico
806	open(21,file=trim(fich_rico),form='formatted')
807	do k=1,llm
808	zlay(k)=0.
809	enddo
810
811	read(21,*) ps_rico,ts_rico
812	prico(1)=ps_rico
813	zrico(1)=0.0
814	do l=2,nlev_rico
815	read(21,*) k,prico(l),zrico(l)
816	enddo
817	close(21)
818
819	do k=1,llm
820	do l=1,80
821	if(prico(l)>play(k)) then
822	if(play(k)>prico(l+1)) then
823	zlay(k)=zrico(l)+(play(k)-prico(l)) * &
824	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
825	else
826	zlay(k)=zrico(l)+(play(k)-prico(80))* &
827	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
828	endif
829	endif
830	enddo
831	print*,k,zlay(k)
832	! U
833	if(0 < zlay(k) .and. zlay(k) < 4000) then
834	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
835	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
836	u_rico(k)= -1.9 + (30.0 + 1.9) / &
837	& (12000 - 4000) * (zlay(k) - 4000)
838	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
839	u_rico(k)=30.0
840	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
841	u_rico(k)=30.0 - (30.0) / &
842	& (20000 - 13000) * (zlay(k) - 13000)
843	else
844	u_rico(k)=0.0
845	endif
846
847	!Q_v
848	if(0 < zlay(k) .and. zlay(k) < 740) then
849	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
850	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
851	q_rico(k)=13.8 + (2.4 - 13.8) / &
852	& (3260 - 740) * (zlay(k) - 740)
853	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
854	q_rico(k)=2.4 + (1.8 - 2.4) / &
855	& (4000 - 3260) * (zlay(k) - 3260)
856	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
857	q_rico(k)=1.8 + (0 - 1.8) / &
858	& (9000 - 4000) * (zlay(k) - 4000)
859	else
860	q_rico(k)=0.0
861	endif
862
863	!T
864	if(0 < zlay(k) .and. zlay(k) < 740) then
865	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
866	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
867	t_rico(k)=292.0 + (278.0 - 292.0) / &
868	& (4000 - 740) * (zlay(k) - 740)
869	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
870	t_rico(k)=278.0 + (203.0 - 278.0) / &
871	& (15000 - 4000) * (zlay(k) - 4000)
872	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
873	t_rico(k)=203.0 + (194.0 - 203.0) / &
874	& (17500 - 15000)* (zlay(k) - 15000)
875	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
876	t_rico(k)=194.0 + (206.0 - 194.0) / &
877	& (20000 - 17500)* (zlay(k) - 17500)
878	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
879	t_rico(k)=206.0 + (270.0 - 206.0) / &
880	& (60000 - 20000)* (zlay(k) - 20000)
881	endif
882
883	! W
884	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
885	w_rico(k)=- (0.005/2260) * zlay(k)
886	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
887	w_rico(k)=- 0.005
888	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
889	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
890	else
891	w_rico(k)=0.0
892	endif
893
894	! dThrz+dTsw0+dTlw0
895	if(0 < zlay(k) .and. zlay(k) < 4000) then
896	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ &
897	& (864004000) zlay(k)
898	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
899	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) / &
900	& (86400(5000 - 4000)) (zlay(k) - 4000)
901	else
902	dth_dyn(k)=0.0
903	endif
904	! dQhrz
905	if(0 < zlay(k) .and. zlay(k) < 3000) then
906	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/ &
907	& (864003000) (zlay(k))
908	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
909	dqh_dyn(k)=0.345 / 86400
910	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
911	dqh_dyn(k)=0.345 / 86400 + &
912	& (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
913	else
914	dqh_dyn(k)=0.0
915	endif
916
917	!? if(play(k)>6e4) then
918	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
919	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
920	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
921	!? *(6e4-play(k))/(6e4-3e4)
922	!? else
923	!? ratqs0(1,k)=ratqshaut
924	!? endif
925
926	enddo
927
928	do k=1,llm
929	q_rico(k)=q_rico(k)/1e3
930	dqh_dyn(k)=dqh_dyn(k)/1e3
931	v_rico(k)=-3.8
932	enddo
933
934	return
935	end
936
937	!======================================================================
938	SUBROUTINE interp_sandu_time(day,day1,annee_ref &
939	& ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu &
940	& ,nlev_sandu,ts_sandu,ts_prof)
941	implicit none
942
943	!---------------------------------------------------------------------------------------
944	! Time interpolation of a 2D field to the timestep corresponding to day
945	!
946	! day: current julian day (e.g. 717538.2)
947	! day1: first day of the simulation
948	! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
949	! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
950	!---------------------------------------------------------------------------------------
951	! inputs:
952	integer annee_ref
953	integer nt_sandu,nlev_sandu
954	integer year_ini_sandu
955	real day, day1,day_ini_sandu,dt_sandu
956	real ts_sandu(nt_sandu)
957	! outputs:
958	real ts_prof
959	! local:
960	integer it_sandu1, it_sandu2
961	real timeit,time_sandu1,time_sandu2,frac
962	! Check that initial day of the simulation consistent with SANDU period:
963	if (annee_ref.ne.2006 ) then
964	print*,'Pour SANDUREF, annee_ref doit etre 2006 '
965	print*,'Changer annee_ref dans run.def'
966	stop
967	endif
968	! if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
969	! print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
970	! print*,'Changer dayref dans run.def'
971	! stop
972	! endif
973
974	! Determine timestep relative to the 1st day of TOGA-COARE:
975	! timeit=(day-day1)*86400.
976	! if (annee_ref.eq.1992) then
977	! timeit=(day-day_ini_sandu)*86400.
978	! else
979	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
980	! endif
981	timeit=(day-day_ini_sandu)*86400
982
983	! Determine the closest observation times:
984	it_sandu1=INT(timeit/dt_sandu)+1
985	it_sandu2=it_sandu1 + 1
986	time_sandu1=(it_sandu1-1)*dt_sandu
987	time_sandu2=(it_sandu2-1)*dt_sandu
988	print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
989	print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2', &
990	& it_sandu1,it_sandu2,time_sandu1,time_sandu2
991
992	if (it_sandu1 .ge. nt_sandu) then
993	write(,) 'PB-stop: day, it_sandu1, it_sandu2, timeit: ' &
994	& ,day,it_sandu1,it_sandu2,timeit/86400.
995	stop
996	endif
997
998	! time interpolation:
999	frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
1000	frac=max(frac,0.0)
1001
1002	ts_prof = ts_sandu(it_sandu2) &
1003	& -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
1004
1005	print*, &
1006	&'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:', &
1007	&day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1, &
1008	&it_sandu2,ts_prof
1009
1010	return
1011	END
1012	!=====================================================================
1013	!-------------------------------------------------------------------------
1014	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu, &
1015	& sens,flat,adv_theta,rad_theta,adv_qt)
1016	implicit none
1017
1018	!-------------------------------------------------------------------------
1019	! Read ARM_CU case forcing data
1020	!-------------------------------------------------------------------------
1021
1022	integer nlev_armcu,nt_armcu
1023	real sens(nt_armcu),flat(nt_armcu)
1024	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
1025	character*80 fich_armcu
1026
1027	integer ip
1028
1029	integer iy,im,id,ih,in
1030
1031	print*,'nlev_armcu',nlev_armcu
1032
1033	open(21,file=trim(fich_armcu),form='formatted')
1034	read(21,'(a)')
1035	do ip = 1, nt_armcu
1036	read(21,'(a)')
1037	read(21,'(a)')
1038	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip), &
1039	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
1040	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip), &
1041	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
1042	enddo
1043	close(21)
1044
1045	223 format(5i3,5f8.3)
1046
1047	return
1048	end
1049
1050	!=====================================================================
1051	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof &
1052	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
1053	& ,ht_prof,vt_prof,hq_prof,vq_prof &
1054	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
1055	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
1056
1057	implicit none
1058
1059	INCLUDE "dimensions.h"
1060
1061	!-------------------------------------------------------------------------
1062	! Vertical interpolation of TOGA-COARE forcing data onto model levels
1063	!-------------------------------------------------------------------------
1064
1065	integer nlevmax
1066	parameter (nlevmax=41)
1067	integer nlev_toga,mxcalc
1068	! real play(llm), plev_prof(nlevmax)
1069	! real t_prof(nlevmax),q_prof(nlevmax)
1070	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1071	! real ht_prof(nlevmax),vt_prof(nlevmax)
1072	! real hq_prof(nlevmax),vq_prof(nlevmax)
1073
1074	real play(llm), plev_prof(nlev_toga)
1075	real t_prof(nlev_toga),q_prof(nlev_toga)
1076	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
1077	real ht_prof(nlev_toga),vt_prof(nlev_toga)
1078	real hq_prof(nlev_toga),vq_prof(nlev_toga)
1079
1080	real t_mod(llm),q_mod(llm)
1081	real u_mod(llm),v_mod(llm), w_mod(llm)
1082	real ht_mod(llm),vt_mod(llm)
1083	real hq_mod(llm),vq_mod(llm)
1084
1085	integer l,k,k1,k2
1086	real frac,frac1,frac2,fact
1087
1088	do l = 1, llm
1089
1090	if (play(l).ge.plev_prof(nlev_toga)) then
1091
1092	mxcalc=l
1093	k1=0
1094	k2=0
1095
1096	if (play(l).le.plev_prof(1)) then
1097
1098	do k = 1, nlev_toga-1
1099	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
1100	k1=k
1101	k2=k+1
1102	endif
1103	enddo
1104
1105	if (k1.eq.0 .or. k2.eq.0) then
1106	write(,) 'PB! k1, k2 = ',k1,k2
1107	write(,) 'l,play(l) = ',l,play(l)/100
1108	do k = 1, nlev_toga-1
1109	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1110	enddo
1111	endif
1112
1113	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1114	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
1115	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
1116	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1117	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1118	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1119	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
1120	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
1121	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
1122	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
1123
1124	else !play>plev_prof(1)
1125
1126	k1=1
1127	k2=2
1128	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1129	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1130	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
1131	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
1132	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
1133	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
1134	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
1135	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
1136	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
1137	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
1138	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
1139
1140	endif ! play.le.plev_prof(1)
1141
1142	else ! above max altitude of forcing file
1143
1144	!jyg
1145	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
1146	fact = max(fact,0.) !jyg
1147	fact = exp(-fact) !jyg
1148	t_mod(l)= t_prof(nlev_toga) !jyg
1149	q_mod(l)= q_prof(nlev_toga)*fact !jyg
1150	u_mod(l)= u_prof(nlev_toga)*fact !jyg
1151	v_mod(l)= v_prof(nlev_toga)*fact !jyg
1152	w_mod(l)= 0.0 !jyg
1153	ht_mod(l)= ht_prof(nlev_toga) !jyg
1154	vt_mod(l)= vt_prof(nlev_toga) !jyg
1155	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
1156	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
1157
1158	endif ! play
1159
1160	enddo ! l
1161
1162	! do l = 1,llm
1163	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
1164	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
1165	! enddo
1166
1167	return
1168	end
1169
1170	!=====================================================================
1171	SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas &
1172	& ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &
1173	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
1174	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
1175	& ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &
1176	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
1177	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
1178
1179	implicit none
1180
1181	INCLUDE "dimensions.h"
1182
1183	!-------------------------------------------------------------------------
1184	! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
1185	!-------------------------------------------------------------------------
1186
1187	integer nlevmax
1188	parameter (nlevmax=41)
1189	integer nlev_cas,mxcalc
1190	! real play(llm), plev_prof(nlevmax)
1191	! real t_prof(nlevmax),q_prof(nlevmax)
1192	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1193	! real ht_prof(nlevmax),vt_prof(nlevmax)
1194	! real hq_prof(nlevmax),vq_prof(nlevmax)
1195
1196	real play(llm), plev_prof_cas(nlev_cas)
1197	real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
1198	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
1199	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
1200	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
1201	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
1202	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
1203	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
1204
1205	real t_mod_cas(llm),q_mod_cas(llm)
1206	real u_mod_cas(llm),v_mod_cas(llm)
1207	real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
1208	real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
1209	real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
1210	real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
1211	real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
1212
1213	integer l,k,k1,k2
1214	real frac,frac1,frac2,fact
1215
1216	do l = 1, llm
1217
1218	if (play(l).ge.plev_prof_cas(nlev_cas)) then
1219
1220	mxcalc=l
1221	k1=0
1222	k2=0
1223
1224	if (play(l).le.plev_prof_cas(1)) then
1225
1226	do k = 1, nlev_cas-1
1227	if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
1228	k1=k
1229	k2=k+1
1230	endif
1231	enddo
1232
1233	if (k1.eq.0 .or. k2.eq.0) then
1234	write(,) 'PB! k1, k2 = ',k1,k2
1235	write(,) 'l,play(l) = ',l,play(l)/100
1236	do k = 1, nlev_cas-1
1237	write(,) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
1238	enddo
1239	endif
1240
1241	frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
1242	t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
1243	q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
1244	u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
1245	v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
1246	ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
1247	vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
1248	w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
1249	du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
1250	hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
1251	vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
1252	dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
1253	hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
1254	vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
1255	dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
1256	ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
1257	vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
1258	dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
1259	hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
1260	vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
1261	dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
1262
1263	else !play>plev_prof_cas(1)
1264
1265	k1=1
1266	k2=2
1267	frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1268	frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1269	t_mod_cas(l)= frac1t_prof_cas(k1) - frac2t_prof_cas(k2)
1270	q_mod_cas(l)= frac1q_prof_cas(k1) - frac2q_prof_cas(k2)
1271	u_mod_cas(l)= frac1u_prof_cas(k1) - frac2u_prof_cas(k2)
1272	v_mod_cas(l)= frac1v_prof_cas(k1) - frac2v_prof_cas(k2)
1273	ug_mod_cas(l)= frac1ug_prof_cas(k1) - frac2ug_prof_cas(k2)
1274	vg_mod_cas(l)= frac1vg_prof_cas(k1) - frac2vg_prof_cas(k2)
1275	w_mod_cas(l)= frac1vitw_prof_cas(k1) - frac2vitw_prof_cas(k2)
1276	du_mod_cas(l)= frac1du_prof_cas(k1) - frac2du_prof_cas(k2)
1277	hu_mod_cas(l)= frac1hu_prof_cas(k1) - frac2hu_prof_cas(k2)
1278	vu_mod_cas(l)= frac1vu_prof_cas(k1) - frac2vu_prof_cas(k2)
1279	dv_mod_cas(l)= frac1dv_prof_cas(k1) - frac2dv_prof_cas(k2)
1280	hv_mod_cas(l)= frac1hv_prof_cas(k1) - frac2hv_prof_cas(k2)
1281	vv_mod_cas(l)= frac1vv_prof_cas(k1) - frac2vv_prof_cas(k2)
1282	dt_mod_cas(l)= frac1dt_prof_cas(k1) - frac2dt_prof_cas(k2)
1283	ht_mod_cas(l)= frac1ht_prof_cas(k1) - frac2ht_prof_cas(k2)
1284	vt_mod_cas(l)= frac1vt_prof_cas(k1) - frac2vt_prof_cas(k2)
1285	dq_mod_cas(l)= frac1dq_prof_cas(k1) - frac2dq_prof_cas(k2)
1286	hq_mod_cas(l)= frac1hq_prof_cas(k1) - frac2hq_prof_cas(k2)
1287	vq_mod_cas(l)= frac1vq_prof_cas(k1) - frac2vq_prof_cas(k2)
1288	dtrad_mod_cas(l)= frac1dtrad_prof_cas(k1) - frac2dtrad_prof_cas(k2)
1289
1290	endif ! play.le.plev_prof_cas(1)
1291
1292	else ! above max altitude of forcing file
1293
1294	!jyg
1295	fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
1296	fact = max(fact,0.) !jyg
1297	fact = exp(-fact) !jyg
1298	t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg
1299	q_mod_cas(l)= q_prof_cas(nlev_cas)*fact !jyg
1300	u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg
1301	v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg
1302	ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg
1303	vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg
1304	w_mod_cas(l)= 0.0 !jyg
1305	du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
1306	hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg
1307	vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg
1308	dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
1309	hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg
1310	vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg
1311	dt_mod_cas(l)= dt_prof_cas(nlev_cas)
1312	ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg
1313	vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg
1314	dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
1315	hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg
1316	vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg
1317	dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg
1318
1319	endif ! play
1320
1321	enddo ! l
1322
1323	! do l = 1,llm
1324	! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
1325	! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
1326	! enddo
1327
1328	return
1329	end
1330	!*****************************************************************************
1331	!=====================================================================
1332	SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof &
1333	& ,th_prof,qv_prof,u_prof,v_prof,o3_prof &
1334	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof &
1335	& ,th_mod,qv_mod,u_mod,v_mod,o3_mod &
1336	& ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
1337
1338	implicit none
1339
1340	INCLUDE "dimensions.h"
1341
1342	!-------------------------------------------------------------------------
1343	! Vertical interpolation of Dice forcing data onto model levels
1344	!-------------------------------------------------------------------------
1345
1346	integer nlevmax
1347	parameter (nlevmax=41)
1348	integer nlev_dice,mxcalc,nt_dice
1349
1350	real play(llm), plev_prof(nlev_dice)
1351	real th_prof(nlev_dice),qv_prof(nlev_dice)
1352	real u_prof(nlev_dice),v_prof(nlev_dice)
1353	real o3_prof(nlev_dice)
1354	real ht_prof(nlev_dice),hq_prof(nlev_dice)
1355	real hu_prof(nlev_dice),hv_prof(nlev_dice)
1356	real w_prof(nlev_dice),omega_prof(nlev_dice)
1357
1358	real th_mod(llm),qv_mod(llm)
1359	real u_mod(llm),v_mod(llm), o3_mod(llm)
1360	real ht_mod(llm),hq_mod(llm)
1361	real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
1362
1363	integer l,k,k1,k2,kp
1364	real aa,frac,frac1,frac2,fact
1365
1366	do l = 1, llm
1367
1368	if (play(l).ge.plev_prof(nlev_dice)) then
1369
1370	mxcalc=l
1371	k1=0
1372	k2=0
1373
1374	if (play(l).le.plev_prof(1)) then
1375
1376	do k = 1, nlev_dice-1
1377	if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
1378	k1=k
1379	k2=k+1
1380	endif
1381	enddo
1382
1383	if (k1.eq.0 .or. k2.eq.0) then
1384	write(,) 'PB! k1, k2 = ',k1,k2
1385	write(,) 'l,play(l) = ',l,play(l)/100
1386	do k = 1, nlev_dice-1
1387	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1388	enddo
1389	endif
1390
1391	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1392	th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
1393	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
1394	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1395	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1396	o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
1397	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
1398	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
1399	hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
1400	hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
1401	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1402	omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
1403
1404	else !play>plev_prof(1)
1405
1406	k1=1
1407	k2=2
1408	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1409	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1410	th_mod(l)= frac1th_prof(k1) - frac2th_prof(k2)
1411	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
1412	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
1413	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
1414	o3_mod(l)= frac1o3_prof(k1) - frac2o3_prof(k2)
1415	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
1416	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
1417	hu_mod(l)= frac1hu_prof(k1) - frac2hu_prof(k2)
1418	hv_mod(l)= frac1hv_prof(k1) - frac2hv_prof(k2)
1419	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
1420	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
1421
1422	endif ! play.le.plev_prof(1)
1423
1424	else ! above max altitude of forcing file
1425
1426	!jyg
1427	fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
1428	fact = max(fact,0.) !jyg
1429	fact = exp(-fact) !jyg
1430	th_mod(l)= th_prof(nlev_dice) !jyg
1431	qv_mod(l)= qv_prof(nlev_dice)*fact !jyg
1432	u_mod(l)= u_prof(nlev_dice)*fact !jyg
1433	v_mod(l)= v_prof(nlev_dice)*fact !jyg
1434	o3_mod(l)= o3_prof(nlev_dice)*fact !jyg
1435	ht_mod(l)= ht_prof(nlev_dice) !jyg
1436	hq_mod(l)= hq_prof(nlev_dice)*fact !jyg
1437	hu_mod(l)= hu_prof(nlev_dice) !jyg
1438	hv_mod(l)= hv_prof(nlev_dice) !jyg
1439	w_mod(l)= 0. !jyg
1440	omega_mod(l)= 0. !jyg
1441
1442	endif ! play
1443
1444	enddo ! l
1445
1446	! do l = 1,llm
1447	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
1448	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
1449	! enddo
1450
1451	return
1452	end
1453
1454	!======================================================================
1455	SUBROUTINE interp_astex_time(day,day1,annee_ref &
1456	& ,year_ini_astex,day_ini_astex,nt_astex,dt_astex &
1457	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
1458	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
1459	& ,ufa_prof,vfa_prof)
1460	implicit none
1461
1462	!---------------------------------------------------------------------------------------
1463	! Time interpolation of a 2D field to the timestep corresponding to day
1464	!
1465	! day: current julian day (e.g. 717538.2)
1466	! day1: first day of the simulation
1467	! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
1468	! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
1469	!---------------------------------------------------------------------------------------
1470
1471	! inputs:
1472	integer annee_ref
1473	integer nt_astex,nlev_astex
1474	integer year_ini_astex
1475	real day, day1,day_ini_astex,dt_astex
1476	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
1477	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
1478	! outputs:
1479	real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
1480	! local:
1481	integer it_astex1, it_astex2
1482	real timeit,time_astex1,time_astex2,frac
1483
1484	! Check that initial day of the simulation consistent with ASTEX period:
1485	if (annee_ref.ne.1992 ) then
1486	print*,'Pour Astex, annee_ref doit etre 1992 '
1487	print*,'Changer annee_ref dans run.def'
1488	stop
1489	endif
1490	if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
1491	print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
1492	print*,'Changer dayref dans run.def'
1493	stop
1494	endif
1495
1496	! Determine timestep relative to the 1st day of TOGA-COARE:
1497	! timeit=(day-day1)*86400.
1498	! if (annee_ref.eq.1992) then
1499	! timeit=(day-day_ini_astex)*86400.
1500	! else
1501	! timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
1502	! endif
1503	timeit=(day-day_ini_astex)*86400
1504
1505	! Determine the closest observation times:
1506	it_astex1=INT(timeit/dt_astex)+1
1507	it_astex2=it_astex1 + 1
1508	time_astex1=(it_astex1-1)*dt_astex
1509	time_astex2=(it_astex2-1)*dt_astex
1510	print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
1511	print *,'it_astex1,it_astex2,time_astex1,time_astex2', &
1512	& it_astex1,it_astex2,time_astex1,time_astex2
1513
1514	if (it_astex1 .ge. nt_astex) then
1515	write(,) 'PB-stop: day, it_astex1, it_astex2, timeit: ' &
1516	& ,day,it_astex1,it_astex2,timeit/86400.
1517	stop
1518	endif
1519
1520	! time interpolation:
1521	frac=(time_astex2-timeit)/(time_astex2-time_astex1)
1522	frac=max(frac,0.0)
1523
1524	div_prof = div_astex(it_astex2) &
1525	& -frac*(div_astex(it_astex2)-div_astex(it_astex1))
1526	ts_prof = ts_astex(it_astex2) &
1527	& -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
1528	ug_prof = ug_astex(it_astex2) &
1529	& -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
1530	vg_prof = vg_astex(it_astex2) &
1531	& -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
1532	ufa_prof = ufa_astex(it_astex2) &
1533	& -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
1534	vfa_prof = vfa_astex(it_astex2) &
1535	& -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
1536
1537	print*, &
1538	&'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:', &
1539	&day,annee_ref,day_ini_astex,timeit/86400.,it_astex1, &
1540	&it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
1541
1542	return
1543	END
1544
1545	!======================================================================
1546	SUBROUTINE interp_toga_time(day,day1,annee_ref &
1547	& ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga &
1548	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
1549	& ,ht_toga,vt_toga,hq_toga,vq_toga &
1550	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
1551	& ,ht_prof,vt_prof,hq_prof,vq_prof)
1552	implicit none
1553
1554	!---------------------------------------------------------------------------------------
1555	! Time interpolation of a 2D field to the timestep corresponding to day
1556	!
1557	! day: current julian day (e.g. 717538.2)
1558	! day1: first day of the simulation
1559	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
1560	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
1561	!---------------------------------------------------------------------------------------
1562
1563	INCLUDE "compar1d.h"
1564
1565	! inputs:
1566	integer annee_ref
1567	integer nt_toga,nlev_toga
1568	integer year_ini_toga
1569	real day, day1,day_ini_toga,dt_toga
1570	real ts_toga(nt_toga)
1571	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
1572	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
1573	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
1574	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1575	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1576	! outputs:
1577	real ts_prof
1578	real plev_prof(nlev_toga),t_prof(nlev_toga)
1579	real q_prof(nlev_toga),u_prof(nlev_toga)
1580	real v_prof(nlev_toga),w_prof(nlev_toga)
1581	real ht_prof(nlev_toga),vt_prof(nlev_toga)
1582	real hq_prof(nlev_toga),vq_prof(nlev_toga)
1583	! local:
1584	integer it_toga1, it_toga2,k
1585	real timeit,time_toga1,time_toga2,frac
1586
1587
1588	if (forcing_type.eq.2) then
1589	! Check that initial day of the simulation consistent with TOGA-COARE period:
1590	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
1591	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
1592	print*,'Changer annee_ref dans run.def'
1593	stop
1594	endif
1595	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
1596	print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'
1597	print*,'Changer dayref dans run.def'
1598	stop
1599	endif
1600	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
1601	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
1602	print*,'Changer dayref ou nday dans run.def'
1603	stop
1604	endif
1605
1606	else if (forcing_type.eq.4) then
1607
1608	! Check that initial day of the simulation consistent with TWP-ICE period:
1609	if (annee_ref.ne.2006) then
1610	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
1611	print*,'Changer annee_ref dans run.def'
1612	stop
1613	endif
1614	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
1615	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
1616	print*,'Changer dayref dans run.def'
1617	stop
1618	endif
1619	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
1620	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
1621	print*,'Changer dayref ou nday dans run.def'
1622	stop
1623	endif
1624
1625	endif
1626
1627	! Determine timestep relative to the 1st day of TOGA-COARE:
1628	! timeit=(day-day1)*86400.
1629	! if (annee_ref.eq.1992) then
1630	! timeit=(day-day_ini_toga)*86400.
1631	! else
1632	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
1633	! endif
1634	timeit=(day-day_ini_toga)*86400
1635
1636	! Determine the closest observation times:
1637	it_toga1=INT(timeit/dt_toga)+1
1638	it_toga2=it_toga1 + 1
1639	time_toga1=(it_toga1-1)*dt_toga
1640	time_toga2=(it_toga2-1)*dt_toga
1641
1642	if (it_toga1 .ge. nt_toga) then
1643	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: ' &
1644	& ,day,it_toga1,it_toga2,timeit/86400.
1645	stop
1646	endif
1647
1648	! time interpolation:
1649	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
1650	frac=max(frac,0.0)
1651
1652	ts_prof = ts_toga(it_toga2) &
1653	& -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
1654
1655	! print*,
1656	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
1657	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
1658
1659	do k=1,nlev_toga
1660	plev_prof(k) = 100.*(plev_toga(k,it_toga2) &
1661	& -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
1662	t_prof(k) = t_toga(k,it_toga2) &
1663	& -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
1664	q_prof(k) = q_toga(k,it_toga2) &
1665	& -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
1666	u_prof(k) = u_toga(k,it_toga2) &
1667	& -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
1668	v_prof(k) = v_toga(k,it_toga2) &
1669	& -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
1670	w_prof(k) = w_toga(k,it_toga2) &
1671	& -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
1672	ht_prof(k) = ht_toga(k,it_toga2) &
1673	& -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
1674	vt_prof(k) = vt_toga(k,it_toga2) &
1675	& -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
1676	hq_prof(k) = hq_toga(k,it_toga2) &
1677	& -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
1678	vq_prof(k) = vq_toga(k,it_toga2) &
1679	& -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
1680	enddo
1681
1682	return
1683	END
1684
1685	!======================================================================
1686	SUBROUTINE interp_dice_time(day,day1,annee_ref &
1687	& ,year_ini_dice,day_ini_dice,nt_dice,dt_dice &
1688	& ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice &
1689	& ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice &
1690	& ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
1691	& ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof &
1692	& ,ustar_prof,psurf_prof,ug_prof,vg_prof &
1693	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
1694	implicit none
1695
1696	!---------------------------------------------------------------------------------------
1697	! Time interpolation of a 2D field to the timestep corresponding to day
1698	!
1699	! day: current julian day (e.g. 717538.2)
1700	! day1: first day of the simulation
1701	! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
1702	! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
1703	!---------------------------------------------------------------------------------------
1704
1705	INCLUDE "compar1d.h"
1706
1707	! inputs:
1708	integer annee_ref
1709	integer nt_dice,nlev_dice
1710	integer year_ini_dice
1711	real day, day1,day_ini_dice,dt_dice
1712	real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
1713	real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
1714	real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
1715	real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
1716	real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
1717	real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
1718	! outputs:
1719	real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
1720	real ustar_prof,psurf_prof,ug_prof,vg_prof
1721	real ht_prof(nlev_dice),hq_prof(nlev_dice)
1722	real hu_prof(nlev_dice),hv_prof(nlev_dice)
1723	real w_prof(nlev_dice),omega_prof(nlev_dice)
1724	! local:
1725	integer it_dice1, it_dice2,k
1726	real timeit,time_dice1,time_dice2,frac
1727
1728
1729	if (forcing_type.eq.7) then
1730	! Check that initial day of the simulation consistent with Dice period:
1731	print *,'annee_ref=',annee_ref
1732	print *,'day1=',day1
1733	print *,'day_ini_dice=',day_ini_dice
1734	if (annee_ref.ne.1999) then
1735	print*,'Pour Dice, annee_ref doit etre 1999'
1736	print*,'Changer annee_ref dans run.def'
1737	stop
1738	endif
1739	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
1740	print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
1741	print*,'Changer dayref dans run.def',day1,day_ini_dice
1742	stop
1743	endif
1744	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
1745	print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
1746	print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
1747	stop
1748	endif
1749
1750	endif
1751
1752	! Determine timestep relative to the 1st day of TOGA-COARE:
1753	! timeit=(day-day1)*86400.
1754	! if (annee_ref.eq.1992) then
1755	! timeit=(day-day_ini_dice)*86400.
1756	! else
1757	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
1758	! endif
1759	timeit=(day-day_ini_dice)*86400
1760
1761	! Determine the closest observation times:
1762	it_dice1=INT(timeit/dt_dice)+1
1763	it_dice2=it_dice1 + 1
1764	time_dice1=(it_dice1-1)*dt_dice
1765	time_dice2=(it_dice2-1)*dt_dice
1766
1767	if (it_dice1 .ge. nt_dice) then
1768	write(,) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
1769	stop
1770	endif
1771
1772	! time interpolation:
1773	frac=(time_dice2-timeit)/(time_dice2-time_dice1)
1774	frac=max(frac,0.0)
1775
1776	shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
1777	lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
1778	lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
1779	swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
1780	tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
1781	ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
1782	psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
1783	ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
1784	vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
1785
1786	! print*,
1787	! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
1788	! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
1789
1790	do k=1,nlev_dice
1791	ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
1792	hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
1793	hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
1794	hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
1795	w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
1796	omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
1797	enddo
1798
1799	return
1800	END
1801
1802	!======================================================================
1803	SUBROUTINE interp_gabls4_time(day,day1,annee_ref &
1804	& ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 &
1805	& ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 &
1806	& ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
1807	implicit none
1808
1809	!---------------------------------------------------------------------------------------
1810	! Time interpolation of a 2D field to the timestep corresponding to day
1811	!
1812	! day: current julian day
1813	! day1: first day of the simulation
1814	! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
1815	! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
1816	!---------------------------------------------------------------------------------------
1817
1818	INCLUDE "compar1d.h"
1819
1820	! inputs:
1821	integer annee_ref
1822	integer nt_gabls4,nlev_gabls4
1823	integer year_ini_gabls4
1824	real day, day1,day_ini_gabls4,dt_gabls4
1825	real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
1826	real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
1827	real tg_gabls4(nt_gabls4), tg_prof
1828	! outputs:
1829	real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
1830	real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
1831	! local:
1832	integer it_gabls41, it_gabls42,k
1833	real timeit,time_gabls41,time_gabls42,frac
1834
1835
1836
1837	! Check that initial day of the simulation consistent with gabls4 period:
1838	if (forcing_type.eq.8 ) then
1839	print *,'annee_ref=',annee_ref
1840	print *,'day1=',day1
1841	print *,'day_ini_gabls4=',day_ini_gabls4
1842	if (annee_ref.ne.2009) then
1843	print*,'Pour gabls4, annee_ref doit etre 2009'
1844	print*,'Changer annee_ref dans run.def'
1845	stop
1846	endif
1847	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
1848	print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
1849	print*,'Changer dayref dans run.def',day1,day_ini_gabls4
1850	stop
1851	endif
1852	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
1853	print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
1854	print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
1855	stop
1856	endif
1857	endif
1858
1859	timeit=(day-day_ini_gabls4)*86400
1860	print *,'day,day_ini_gabls4=',day,day_ini_gabls4
1861	print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
1862
1863	! Determine the closest observation times:
1864	it_gabls41=INT(timeit/dt_gabls4)+1
1865	it_gabls42=it_gabls41 + 1
1866	time_gabls41=(it_gabls41-1)*dt_gabls4
1867	time_gabls42=(it_gabls42-1)*dt_gabls4
1868
1869	if (it_gabls41 .ge. nt_gabls4) then
1870	write(,) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
1871	stop
1872	endif
1873
1874	! time interpolation:
1875	frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
1876	frac=max(frac,0.0)
1877
1878
1879	do k=1,nlev_gabls4
1880	ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
1881	vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
1882	ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
1883	hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
1884	enddo
1885	tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
1886	return
1887	END
1888
1889	!======================================================================
1890	SUBROUTINE interp_armcu_time(day,day1,annee_ref &
1891	& ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu &
1892	& ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu &
1893	& ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
1894	implicit none
1895
1896	!---------------------------------------------------------------------------------------
1897	! Time interpolation of a 2D field to the timestep corresponding to day
1898	!
1899	! day: current julian day (e.g. 717538.2)
1900	! day1: first day of the simulation
1901	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
1902	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
1903	! fs= sensible flux
1904	! fl= latent flux
1905	! at,rt,aqt= advective and radiative tendencies
1906	!---------------------------------------------------------------------------------------
1907
1908	! inputs:
1909	integer annee_ref
1910	integer nt_armcu,nlev_armcu
1911	integer year_ini_armcu
1912	real day, day1,day_ini_armcu,dt_armcu
1913	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
1914	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
1915	! outputs:
1916	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
1917	! local:
1918	integer it_armcu1, it_armcu2,k
1919	real timeit,time_armcu1,time_armcu2,frac
1920
1921	! Check that initial day of the simulation consistent with ARMCU period:
1922	if (annee_ref.ne.1997 ) then
1923	print*,'Pour ARMCU, annee_ref doit etre 1997 '
1924	print*,'Changer annee_ref dans run.def'
1925	stop
1926	endif
1927
1928	timeit=(day-day_ini_armcu)*86400
1929
1930	! Determine the closest observation times:
1931	it_armcu1=INT(timeit/dt_armcu)+1
1932	it_armcu2=it_armcu1 + 1
1933	time_armcu1=(it_armcu1-1)*dt_armcu
1934	time_armcu2=(it_armcu2-1)*dt_armcu
1935	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
1936	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', &
1937	& it_armcu1,it_armcu2,time_armcu1,time_armcu2
1938
1939	if (it_armcu1 .ge. nt_armcu) then
1940	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' &
1941	& ,day,it_armcu1,it_armcu2,timeit/86400.
1942	stop
1943	endif
1944
1945	! time interpolation:
1946	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
1947	frac=max(frac,0.0)
1948
1949	fs_prof = fs_armcu(it_armcu2) &
1950	& -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
1951	fl_prof = fl_armcu(it_armcu2) &
1952	& -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
1953	at_prof = at_armcu(it_armcu2) &
1954	& -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
1955	rt_prof = rt_armcu(it_armcu2) &
1956	& -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
1957	aqt_prof = aqt_armcu(it_armcu2) &
1958	& -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
1959
1960	print*, &
1961	&'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', &
1962	&day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, &
1963	&it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
1964
1965	return
1966	END
1967
1968	!=====================================================================
1969	subroutine readprofiles(nlev_max,kmax,ntrac,height, &
1970	& thlprof,qtprof,uprof, &
1971	& vprof,e12prof,ugprof,vgprof, &
1972	& wfls,dqtdxls,dqtdyls,dqtdtls, &
1973	& thlpcar,tracer,nt1,nt2)
1974	implicit none
1975
1976	integer nlev_max,kmax,kmax2,ntrac
1977	logical :: llesread = .true.
1978
1979	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), &
1980	& uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), &
1981	& ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), &
1982	& dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), &
1983	& thlpcar(nlev_max),tracer(nlev_max,ntrac)
1984
1985	real height1(nlev_max)
1986
1987	integer, parameter :: ilesfile=1
1988	integer :: ierr,k,itrac,nt1,nt2
1989
1990	if(.not.(llesread)) return
1991
1992	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
1993	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
1994	read (ilesfile,*) kmax
1995	do k=1,kmax
1996	read (ilesfile,*) height1(k),thlprof(k),qtprof (k), &
1997	& uprof (k),vprof (k),e12prof(k)
1998	enddo
1999	close(ilesfile)
2000
2001	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2002	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2003	read (ilesfile,*) kmax2
2004	if (kmax .ne. kmax2) then
2005	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2006	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2007	stop 'lecture profiles'
2008	endif
2009	do k=1,kmax
2010	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &
2011	& dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2012	end do
2013	do k=1,kmax
2014	if (height(k) .ne. height1(k)) then
2015	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2016	print *, 'les niveaux different : ',k,height1(k), height(k)
2017	stop
2018	endif
2019	end do
2020	close(ilesfile)
2021
2022	open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
2023	if (ierr /= 0) then
2024	print*,'WARNING : trac.inp does not exist'
2025	else
2026	read (ilesfile,*) kmax2,nt1,nt2
2027	if (nt2>ntrac) then
2028	stop 'Augmenter le nombre de traceurs dans traceur.def'
2029	endif
2030	if (kmax .ne. kmax2) then
2031	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2032	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2033	stop 'lecture profiles'
2034	endif
2035	do k=1,kmax
2036	read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
2037	end do
2038	close(ilesfile)
2039	endif
2040
2041	return
2042	end
2043	!======================================================================
2044	subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &
2045	& thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
2046	!======================================================================
2047	implicit none
2048
2049	integer nlev_max,kmax
2050	logical :: llesread = .true.
2051
2052	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2053	real thlprof(nlev_max)
2054	real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2055	real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
2056
2057	integer, parameter :: ilesfile=1
2058	integer :: k,ierr
2059
2060	if(.not.(llesread)) return
2061
2062	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2063	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2064	read (ilesfile,*) kmax
2065	do k=1,kmax
2066	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
2067	& qprof (k),uprof(k), vprof(k), wprof(k), &
2068	& omega (k),o3mmr(k)
2069	enddo
2070	close(ilesfile)
2071
2072	return
2073	end
2074
2075	!======================================================================
2076	subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &
2077	& thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
2078	!======================================================================
2079	implicit none
2080
2081	integer nlev_max,kmax
2082	logical :: llesread = .true.
2083
2084	real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &
2085	& thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &
2086	& qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &
2087	& wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
2088
2089	integer, parameter :: ilesfile=1
2090	integer :: ierr,k
2091
2092	if(.not.(llesread)) return
2093
2094	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2095	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2096	read (ilesfile,*) kmax
2097	do k=1,kmax
2098	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
2099	& qvprof (k),qlprof (k),qtprof (k), &
2100	& uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
2101	enddo
2102	close(ilesfile)
2103
2104	return
2105	end
2106
2107
2108
2109	!======================================================================
2110	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &
2111	& vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2112	!======================================================================
2113	implicit none
2114
2115	integer nlev_max,kmax
2116	logical :: llesread = .true.
2117
2118	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2119	real thetaprof(nlev_max),rvprof(nlev_max)
2120	real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2121	real aprof(nlev_max+1),bprof(nlev_max+1)
2122
2123	integer, parameter :: ilesfile=1
2124	integer, parameter :: ifile=2
2125	integer :: ierr,jtot,k
2126
2127	if(.not.(llesread)) return
2128
2129	! Read profiles at full levels
2130	IF(nlev_max.EQ.19) THEN
2131	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2132	print *,'On ouvre prof.inp.19'
2133	ELSE
2134	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2135	print *,'On ouvre prof.inp.40'
2136	ENDIF
2137	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2138	read (ilesfile,*) kmax
2139	do k=1,kmax
2140	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &
2141	& thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2142	enddo
2143	close(ilesfile)
2144
2145	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
2146	IF(nlev_max.EQ.19) THEN
2147	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2148	print *,'On ouvre proh.inp.19'
2149	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2150	ELSE
2151	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2152	print *,'On ouvre proh.inp.40'
2153	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2154	ENDIF
2155	read (ifile,*) kmax
2156	do k=1,kmax
2157	read (ifile,*) jtot,aprof(k),bprof(k)
2158	enddo
2159	close(ifile)
2160
2161	return
2162	end
2163
2164	!=====================================================================
2165	subroutine read_fire(fich_fire,nlevel,ntime &
2166	& ,zz,thl,qt,u,v,tke &
2167	& ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
2168
2169	!program reading forcings of the FIRE case study
2170
2171
2172	use netcdf, only: nf90_get_var
2173	implicit none
2174
2175	INCLUDE "netcdf.inc"
2176
2177	integer ntime,nlevel
2178	character*80 :: fich_fire
2179	real*8 zz(nlevel)
2180
2181	real*8 thl(nlevel)
2182	real*8 qt(nlevel),u(nlevel)
2183	real*8 v(nlevel),tke(nlevel)
2184	real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
2185	real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
2186	real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
2187
2188	integer nid, ierr
2189	integer nbvar3d
2190	parameter(nbvar3d=30)
2191	integer var3didin(nbvar3d)
2192
2193	ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
2194	if (ierr.NE.NF_NOERR) then
2195	write(,) 'ERROR: Pb opening forcings nc file '
2196	write(,) NF_STRERROR(ierr)
2197	stop ""
2198	endif
2199
2200
2201	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
2202	if(ierr/=NF_NOERR) then
2203	write(,) NF_STRERROR(ierr)
2204	stop 'lev'
2205	endif
2206
2207
2208	ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
2209	if(ierr/=NF_NOERR) then
2210	write(,) NF_STRERROR(ierr)
2211	stop 'temp'
2212	endif
2213
2214	ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
2215	if(ierr/=NF_NOERR) then
2216	write(,) NF_STRERROR(ierr)
2217	stop 'qv'
2218	endif
2219
2220	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
2221	if(ierr/=NF_NOERR) then
2222	write(,) NF_STRERROR(ierr)
2223	stop 'u'
2224	endif
2225
2226	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
2227	if(ierr/=NF_NOERR) then
2228	write(,) NF_STRERROR(ierr)
2229	stop 'v'
2230	endif
2231
2232	ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
2233	if(ierr/=NF_NOERR) then
2234	write(,) NF_STRERROR(ierr)
2235	stop 'tke'
2236	endif
2237
2238	ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
2239	if(ierr/=NF_NOERR) then
2240	write(,) NF_STRERROR(ierr)
2241	stop 'ug'
2242	endif
2243
2244	ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
2245	if(ierr/=NF_NOERR) then
2246	write(,) NF_STRERROR(ierr)
2247	stop 'vg'
2248	endif
2249
2250	ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
2251	if(ierr/=NF_NOERR) then
2252	write(,) NF_STRERROR(ierr)
2253	stop 'wls'
2254	endif
2255
2256	ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
2257	if(ierr/=NF_NOERR) then
2258	write(,) NF_STRERROR(ierr)
2259	stop 'dqtdx'
2260	endif
2261
2262	ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
2263	if(ierr/=NF_NOERR) then
2264	write(,) NF_STRERROR(ierr)
2265	stop 'dqtdy'
2266	endif
2267
2268	ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
2269	if(ierr/=NF_NOERR) then
2270	write(,) NF_STRERROR(ierr)
2271	stop 'dqtdt'
2272	endif
2273
2274	ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
2275	if(ierr/=NF_NOERR) then
2276	write(,) NF_STRERROR(ierr)
2277	stop 'thl_rad'
2278	endif
2279	!dimensions lecture
2280	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
2281
2282	ierr = NF90_GET_VAR(nid,var3didin(1),zz)
2283	if(ierr/=NF_NOERR) then
2284	write(,) NF_STRERROR(ierr)
2285	stop "getvarup"
2286	endif
2287	! write(,)'lecture z ok',zz
2288
2289	ierr = NF90_GET_VAR(nid,var3didin(2),thl)
2290	if(ierr/=NF_NOERR) then
2291	write(,) NF_STRERROR(ierr)
2292	stop "getvarup"
2293	endif
2294	! write(,)'lecture thl ok',thl
2295
2296	ierr = NF90_GET_VAR(nid,var3didin(3),qt)
2297	if(ierr/=NF_NOERR) then
2298	write(,) NF_STRERROR(ierr)
2299	stop "getvarup"
2300	endif
2301	! write(,)'lecture qt ok',qt
2302
2303	ierr = NF90_GET_VAR(nid,var3didin(4),u)
2304	if(ierr/=NF_NOERR) then
2305	write(,) NF_STRERROR(ierr)
2306	stop "getvarup"
2307	endif
2308	! write(,)'lecture u ok',u
2309
2310	ierr = NF90_GET_VAR(nid,var3didin(5),v)
2311	if(ierr/=NF_NOERR) then
2312	write(,) NF_STRERROR(ierr)
2313	stop "getvarup"
2314	endif
2315	! write(,)'lecture v ok',v
2316
2317	ierr = NF90_GET_VAR(nid,var3didin(6),tke)
2318	if(ierr/=NF_NOERR) then
2319	write(,) NF_STRERROR(ierr)
2320	stop "getvarup"
2321	endif
2322	! write(,)'lecture tke ok',tke
2323
2324	ierr = NF90_GET_VAR(nid,var3didin(7),ug)
2325	if(ierr/=NF_NOERR) then
2326	write(,) NF_STRERROR(ierr)
2327	stop "getvarup"
2328	endif
2329	! write(,)'lecture ug ok',ug
2330
2331	ierr = NF90_GET_VAR(nid,var3didin(8),vg)
2332	if(ierr/=NF_NOERR) then
2333	write(,) NF_STRERROR(ierr)
2334	stop "getvarup"
2335	endif
2336	! write(,)'lecture vg ok',vg
2337
2338	ierr = NF90_GET_VAR(nid,var3didin(9),wls)
2339	if(ierr/=NF_NOERR) then
2340	write(,) NF_STRERROR(ierr)
2341	stop "getvarup"
2342	endif
2343	! write(,)'lecture wls ok',wls
2344
2345	ierr = NF90_GET_VAR(nid,var3didin(10),dqtdx)
2346	if(ierr/=NF_NOERR) then
2347	write(,) NF_STRERROR(ierr)
2348	stop "getvarup"
2349	endif
2350	! write(,)'lecture dqtdx ok',dqtdx
2351
2352	ierr = NF90_GET_VAR(nid,var3didin(11),dqtdy)
2353	if(ierr/=NF_NOERR) then
2354	write(,) NF_STRERROR(ierr)
2355	stop "getvarup"
2356	endif
2357	! write(,)'lecture dqtdy ok',dqtdy
2358
2359	ierr = NF90_GET_VAR(nid,var3didin(12),dqtdt)
2360	if(ierr/=NF_NOERR) then
2361	write(,) NF_STRERROR(ierr)
2362	stop "getvarup"
2363	endif
2364	! write(,)'lecture dqtdt ok',dqtdt
2365
2366	ierr = NF90_GET_VAR(nid,var3didin(13),thl_rad)
2367	if(ierr/=NF_NOERR) then
2368	write(,) NF_STRERROR(ierr)
2369	stop "getvarup"
2370	endif
2371	! write(,)'lecture thl_rad ok',thl_rad
2372
2373	return
2374	end subroutine read_fire
2375	!=====================================================================
2376	subroutine read_dice(fich_dice,nlevel,ntime &
2377	& ,zz,pres,t,qv,u,v,o3 &
2378	& ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg &
2379	& ,hadvt,hadvq,hadvu,hadvv,w,omega)
2380
2381	!program reading initial profils and forcings of the Dice case study
2382
2383	use netcdf, only: nf90_get_var
2384
2385	implicit none
2386
2387	INCLUDE "netcdf.inc"
2388	INCLUDE "YOMCST.h"
2389
2390	integer ntime,nlevel
2391	integer l,k
2392	character*80 :: fich_dice
2393	real*8 time(ntime)
2394	real*8 zz(nlevel)
2395
2396	real*8 th(nlevel),pres(nlevel),t(nlevel)
2397	real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
2398	real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
2399	real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
2400	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
2401	real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
2402	real*8 pzero
2403
2404	integer nid, ierr
2405	integer nbvar3d
2406	parameter(nbvar3d=30)
2407	integer var3didin(nbvar3d)
2408
2409	pzero=100000.
2410	ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
2411	if (ierr.NE.NF_NOERR) then
2412	write(,) 'ERROR: Pb opening forcings nc file '
2413	write(,) NF_STRERROR(ierr)
2414	stop ""
2415	endif
2416
2417
2418	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
2419	if(ierr/=NF_NOERR) then
2420	write(,) NF_STRERROR(ierr)
2421	stop 'height'
2422	endif
2423
2424	ierr=NF_INQ_VARID(nid,"pf",var3didin(11))
2425	if(ierr/=NF_NOERR) then
2426	write(,) NF_STRERROR(ierr)
2427	stop 'pf'
2428	endif
2429
2430	ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
2431	if(ierr/=NF_NOERR) then
2432	write(,) NF_STRERROR(ierr)
2433	stop 'theta'
2434	endif
2435
2436	ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
2437	if(ierr/=NF_NOERR) then
2438	write(,) NF_STRERROR(ierr)
2439	stop 'qv'
2440	endif
2441
2442	ierr=NF_INQ_VARID(nid,"u",var3didin(14))
2443	if(ierr/=NF_NOERR) then
2444	write(,) NF_STRERROR(ierr)
2445	stop 'u'
2446	endif
2447
2448	ierr=NF_INQ_VARID(nid,"v",var3didin(15))
2449	if(ierr/=NF_NOERR) then
2450	write(,) NF_STRERROR(ierr)
2451	stop 'v'
2452	endif
2453
2454	ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
2455	if(ierr/=NF_NOERR) then
2456	write(,) NF_STRERROR(ierr)
2457	stop 'o3'
2458	endif
2459
2460	ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
2461	if(ierr/=NF_NOERR) then
2462	write(,) NF_STRERROR(ierr)
2463	stop 'shf'
2464	endif
2465
2466	ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
2467	if(ierr/=NF_NOERR) then
2468	write(,) NF_STRERROR(ierr)
2469	stop 'lhf'
2470	endif
2471
2472	ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
2473	if(ierr/=NF_NOERR) then
2474	write(,) NF_STRERROR(ierr)
2475	stop 'lwup'
2476	endif
2477
2478	ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
2479	if(ierr/=NF_NOERR) then
2480	write(,) NF_STRERROR(ierr)
2481	stop 'dqtdx'
2482	endif
2483
2484	ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
2485	if(ierr/=NF_NOERR) then
2486	write(,) NF_STRERROR(ierr)
2487	stop 'Tg'
2488	endif
2489
2490	ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
2491	if(ierr/=NF_NOERR) then
2492	write(,) NF_STRERROR(ierr)
2493	stop 'ustar'
2494	endif
2495
2496	ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
2497	if(ierr/=NF_NOERR) then
2498	write(,) NF_STRERROR(ierr)
2499	stop 'psurf'
2500	endif
2501
2502	ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
2503	if(ierr/=NF_NOERR) then
2504	write(,) NF_STRERROR(ierr)
2505	stop 'Ug'
2506	endif
2507
2508	ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
2509	if(ierr/=NF_NOERR) then
2510	write(,) NF_STRERROR(ierr)
2511	stop 'Vg'
2512	endif
2513
2514	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
2515	if(ierr/=NF_NOERR) then
2516	write(,) NF_STRERROR(ierr)
2517	stop 'hadvT'
2518	endif
2519
2520	ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
2521	if(ierr/=NF_NOERR) then
2522	write(,) NF_STRERROR(ierr)
2523	stop 'hadvq'
2524	endif
2525
2526	ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
2527	if(ierr/=NF_NOERR) then
2528	write(,) NF_STRERROR(ierr)
2529	stop 'hadvu'
2530	endif
2531
2532	ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
2533	if(ierr/=NF_NOERR) then
2534	write(,) NF_STRERROR(ierr)
2535	stop 'hadvv'
2536	endif
2537
2538	ierr=NF_INQ_VARID(nid,"w",var3didin(21))
2539	if(ierr/=NF_NOERR) then
2540	write(,) NF_STRERROR(ierr)
2541	stop 'w'
2542	endif
2543
2544	ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
2545	if(ierr/=NF_NOERR) then
2546	write(,) NF_STRERROR(ierr)
2547	stop 'omega'
2548	endif
2549	!dimensions lecture
2550	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
2551
2552	ierr = NF90_GET_VAR(nid,var3didin(1),zz)
2553	if(ierr/=NF_NOERR) then
2554	write(,) NF_STRERROR(ierr)
2555	stop "getvarup"
2556	endif
2557	! write(,)'lecture zz ok',zz
2558
2559	ierr = NF90_GET_VAR(nid,var3didin(11),pres)
2560	if(ierr/=NF_NOERR) then
2561	write(,) NF_STRERROR(ierr)
2562	stop "getvarup"
2563	endif
2564	! write(,)'lecture pres ok',pres
2565
2566	ierr = NF90_GET_VAR(nid,var3didin(12),th)
2567	if(ierr/=NF_NOERR) then
2568	write(,) NF_STRERROR(ierr)
2569	stop "getvarup"
2570	endif
2571	! write(,)'lecture th ok',th
2572	do k=1,nlevel
2573	t(k)=th(k)(pres(k)/pzero)*rkappa
2574	enddo
2575
2576	ierr = NF90_GET_VAR(nid,var3didin(13),qv)
2577	if(ierr/=NF_NOERR) then
2578	write(,) NF_STRERROR(ierr)
2579	stop "getvarup"
2580	endif
2581	! write(,)'lecture qv ok',qv
2582
2583	ierr = NF90_GET_VAR(nid,var3didin(14),u)
2584	if(ierr/=NF_NOERR) then
2585	write(,) NF_STRERROR(ierr)
2586	stop "getvarup"
2587	endif
2588	! write(,)'lecture u ok',u
2589
2590	ierr = NF90_GET_VAR(nid,var3didin(15),v)
2591	if(ierr/=NF_NOERR) then
2592	write(,) NF_STRERROR(ierr)
2593	stop "getvarup"
2594	endif
2595	! write(,)'lecture v ok',v
2596
2597	ierr = NF90_GET_VAR(nid,var3didin(16),o3)
2598	if(ierr/=NF_NOERR) then
2599	write(,) NF_STRERROR(ierr)
2600	stop "getvarup"
2601	endif
2602	! write(,)'lecture o3 ok',o3
2603
2604	ierr = NF90_GET_VAR(nid,var3didin(2),shf)
2605	if(ierr/=NF_NOERR) then
2606	write(,) NF_STRERROR(ierr)
2607	stop "getvarup"
2608	endif
2609	! write(,)'lecture shf ok',shf
2610
2611	ierr = NF90_GET_VAR(nid,var3didin(3),lhf)
2612	if(ierr/=NF_NOERR) then
2613	write(,) NF_STRERROR(ierr)
2614	stop "getvarup"
2615	endif
2616	! write(,)'lecture lhf ok',lhf
2617
2618	ierr = NF90_GET_VAR(nid,var3didin(4),lwup)
2619	if(ierr/=NF_NOERR) then
2620	write(,) NF_STRERROR(ierr)
2621	stop "getvarup"
2622	endif
2623	! write(,)'lecture lwup ok',lwup
2624
2625	ierr = NF90_GET_VAR(nid,var3didin(5),swup)
2626	if(ierr/=NF_NOERR) then
2627	write(,) NF_STRERROR(ierr)
2628	stop "getvarup"
2629	endif
2630	! write(,)'lecture swup ok',swup
2631
2632	ierr = NF90_GET_VAR(nid,var3didin(6),tg)
2633	if(ierr/=NF_NOERR) then
2634	write(,) NF_STRERROR(ierr)
2635	stop "getvarup"
2636	endif
2637	! write(,)'lecture tg ok',tg
2638
2639	ierr = NF90_GET_VAR(nid,var3didin(7),ustar)
2640	if(ierr/=NF_NOERR) then
2641	write(,) NF_STRERROR(ierr)
2642	stop "getvarup"
2643	endif
2644	! write(,)'lecture ustar ok',ustar
2645
2646	ierr = NF90_GET_VAR(nid,var3didin(8),psurf)
2647	if(ierr/=NF_NOERR) then
2648	write(,) NF_STRERROR(ierr)
2649	stop "getvarup"
2650	endif
2651	! write(,)'lecture psurf ok',psurf
2652
2653	ierr = NF90_GET_VAR(nid,var3didin(9),ug)
2654	if(ierr/=NF_NOERR) then
2655	write(,) NF_STRERROR(ierr)
2656	stop "getvarup"
2657	endif
2658	! write(,)'lecture ug ok',ug
2659
2660	ierr = NF90_GET_VAR(nid,var3didin(10),vg)
2661	if(ierr/=NF_NOERR) then
2662	write(,) NF_STRERROR(ierr)
2663	stop "getvarup"
2664	endif
2665	! write(,)'lecture vg ok',vg
2666
2667	ierr = NF90_GET_VAR(nid,var3didin(17),hadvt)
2668	if(ierr/=NF_NOERR) then
2669	write(,) NF_STRERROR(ierr)
2670	stop "getvarup"
2671	endif
2672	! write(,)'lecture hadvt ok',hadvt
2673
2674	ierr = NF90_GET_VAR(nid,var3didin(18),hadvq)
2675	if(ierr/=NF_NOERR) then
2676	write(,) NF_STRERROR(ierr)
2677	stop "getvarup"
2678	endif
2679	! write(,)'lecture hadvq ok',hadvq
2680
2681	ierr = NF90_GET_VAR(nid,var3didin(19),hadvu)
2682	if(ierr/=NF_NOERR) then
2683	write(,) NF_STRERROR(ierr)
2684	stop "getvarup"
2685	endif
2686	! write(,)'lecture hadvu ok',hadvu
2687
2688	ierr = NF90_GET_VAR(nid,var3didin(20),hadvv)
2689	if(ierr/=NF_NOERR) then
2690	write(,) NF_STRERROR(ierr)
2691	stop "getvarup"
2692	endif
2693	! write(,)'lecture hadvv ok',hadvv
2694
2695	ierr = NF90_GET_VAR(nid,var3didin(21),w)
2696	if(ierr/=NF_NOERR) then
2697	write(,) NF_STRERROR(ierr)
2698	stop "getvarup"
2699	endif
2700	! write(,)'lecture w ok',w
2701
2702	ierr = NF90_GET_VAR(nid,var3didin(22),omega)
2703	if(ierr/=NF_NOERR) then
2704	write(,) NF_STRERROR(ierr)
2705	stop "getvarup"
2706	endif
2707	! write(,)'lecture omega ok',omega
2708
2709	return
2710	end subroutine read_dice
2711	!=====================================================================
2712	subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol &
2713	& ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
2714
2715	!program reading initial profils and forcings of the Gabls4 case study
2716
2717	use netcdf, only: nf90_get_var
2718
2719	implicit none
2720
2721	INCLUDE "netcdf.inc"
2722
2723	integer ntime,nlevel,nsol
2724	integer l,k
2725	character*80 :: fich_gabls4
2726	real*8 time(ntime)
2727
2728	! ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees
2729	! dans un ordre inverse par rapport a la convention LMDZ
2730	! ==> il faut tout inverser (MPL 20141024)
2731	! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
2732	real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
2733	real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
2734	real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
2735
2736	real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
2737	real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
2738	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
2739
2740	real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
2741	real*8 tg(ntime)
2742	integer nid, ierr
2743	integer nbvar3d
2744	parameter(nbvar3d=30)
2745	integer var3didin(nbvar3d)
2746
2747	ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
2748	if (ierr.NE.NF_NOERR) then
2749	write(,) 'ERROR: Pb opening forcings nc file '
2750	write(,) NF_STRERROR(ierr)
2751	stop ""
2752	endif
2753
2754
2755	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
2756	if(ierr/=NF_NOERR) then
2757	write(,) NF_STRERROR(ierr)
2758	stop 'height'
2759	endif
2760
2761	ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
2762	if(ierr/=NF_NOERR) then
2763	write(,) NF_STRERROR(ierr)
2764	stop 'depth_sn'
2765	endif
2766
2767	ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
2768	if(ierr/=NF_NOERR) then
2769	write(,) NF_STRERROR(ierr)
2770	stop 'Ug'
2771	endif
2772
2773	ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
2774	if(ierr/=NF_NOERR) then
2775	write(,) NF_STRERROR(ierr)
2776	stop 'Vg'
2777	endif
2778	ierr=NF_INQ_VARID(nid,"pf",var3didin(5))
2779	if(ierr/=NF_NOERR) then
2780	write(,) NF_STRERROR(ierr)
2781	stop 'pf'
2782	endif
2783
2784	ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
2785	if(ierr/=NF_NOERR) then
2786	write(,) NF_STRERROR(ierr)
2787	stop 'theta'
2788	endif
2789
2790	ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
2791	if(ierr/=NF_NOERR) then
2792	write(,) NF_STRERROR(ierr)
2793	stop 'tempe'
2794	endif
2795
2796	ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
2797	if(ierr/=NF_NOERR) then
2798	write(,) NF_STRERROR(ierr)
2799	stop 'qv'
2800	endif
2801
2802	ierr=NF_INQ_VARID(nid,"u",var3didin(9))
2803	if(ierr/=NF_NOERR) then
2804	write(,) NF_STRERROR(ierr)
2805	stop 'u'
2806	endif
2807
2808	ierr=NF_INQ_VARID(nid,"v",var3didin(10))
2809	if(ierr/=NF_NOERR) then
2810	write(,) NF_STRERROR(ierr)
2811	stop 'v'
2812	endif
2813
2814	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
2815	if(ierr/=NF_NOERR) then
2816	write(,) NF_STRERROR(ierr)
2817	stop 'hadvt'
2818	endif
2819
2820	ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
2821	if(ierr/=NF_NOERR) then
2822	write(,) NF_STRERROR(ierr)
2823	stop 'hadvq'
2824	endif
2825
2826	ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
2827	if(ierr/=NF_NOERR) then
2828	write(,) NF_STRERROR(ierr)
2829	stop 'tsnow'
2830	endif
2831
2832	ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
2833	if(ierr/=NF_NOERR) then
2834	write(,) NF_STRERROR(ierr)
2835	stop 'snow_density'
2836	endif
2837
2838	ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
2839	if(ierr/=NF_NOERR) then
2840	write(,) NF_STRERROR(ierr)
2841	stop 'Tg'
2842	endif
2843
2844
2845	!dimensions lecture
2846	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
2847
2848	ierr = NF90_GET_VAR(nid,var3didin(1),zz_i)
2849	if(ierr/=NF_NOERR) then
2850	write(,) NF_STRERROR(ierr)
2851	stop "getvarup"
2852	endif
2853
2854	ierr = NF90_GET_VAR(nid,var3didin(2),depth_sn)
2855	if(ierr/=NF_NOERR) then
2856	write(,) NF_STRERROR(ierr)
2857	stop "getvarup"
2858	endif
2859
2860	ierr = NF90_GET_VAR(nid,var3didin(3),ug_i)
2861	if(ierr/=NF_NOERR) then
2862	write(,) NF_STRERROR(ierr)
2863	stop "getvarup"
2864	endif
2865
2866	ierr = NF90_GET_VAR(nid,var3didin(4),vg_i)
2867	if(ierr/=NF_NOERR) then
2868	write(,) NF_STRERROR(ierr)
2869	stop "getvarup"
2870	endif
2871
2872	ierr = NF90_GET_VAR(nid,var3didin(5),pf_i)
2873	if(ierr/=NF_NOERR) then
2874	write(,) NF_STRERROR(ierr)
2875	stop "getvarup"
2876	endif
2877
2878	ierr = NF90_GET_VAR(nid,var3didin(6),th_i)
2879	if(ierr/=NF_NOERR) then
2880	write(,) NF_STRERROR(ierr)
2881	stop "getvarup"
2882	endif
2883
2884	ierr = NF90_GET_VAR(nid,var3didin(7),t_i)
2885	if(ierr/=NF_NOERR) then
2886	write(,) NF_STRERROR(ierr)
2887	stop "getvarup"
2888	endif
2889
2890	ierr = NF90_GET_VAR(nid,var3didin(8),qv_i)
2891	if(ierr/=NF_NOERR) then
2892	write(,) NF_STRERROR(ierr)
2893	stop "getvarup"
2894	endif
2895
2896	ierr = NF90_GET_VAR(nid,var3didin(9),u_i)
2897	if(ierr/=NF_NOERR) then
2898	write(,) NF_STRERROR(ierr)
2899	stop "getvarup"
2900	endif
2901
2902	ierr = NF90_GET_VAR(nid,var3didin(10),v_i)
2903	if(ierr/=NF_NOERR) then
2904	write(,) NF_STRERROR(ierr)
2905	stop "getvarup"
2906	endif
2907
2908	ierr = NF90_GET_VAR(nid,var3didin(11),hadvt_i)
2909	if(ierr/=NF_NOERR) then
2910	write(,) NF_STRERROR(ierr)
2911	stop "getvarup"
2912	endif
2913
2914	ierr = NF90_GET_VAR(nid,var3didin(12),hadvq_i)
2915	if(ierr/=NF_NOERR) then
2916	write(,) NF_STRERROR(ierr)
2917	stop "getvarup"
2918	endif
2919
2920	ierr = NF90_GET_VAR(nid,var3didin(14),tsnow)
2921	if(ierr/=NF_NOERR) then
2922	write(,) NF_STRERROR(ierr)
2923	stop "getvarup"
2924	endif
2925
2926	ierr = NF90_GET_VAR(nid,var3didin(15),snow_dens)
2927	if(ierr/=NF_NOERR) then
2928	write(,) NF_STRERROR(ierr)
2929	stop "getvarup"
2930	endif
2931
2932	ierr = NF90_GET_VAR(nid,var3didin(16),tg)
2933	if(ierr/=NF_NOERR) then
2934	write(,) NF_STRERROR(ierr)
2935	stop "getvarup"
2936	endif
2937
2938	! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
2939	do k=1,nlevel
2940	zz(k)=zz_i(nlevel+1-k)
2941	ug(k,:)=ug_i(nlevel+1-k,:)
2942	vg(k,:)=vg_i(nlevel+1-k,:)
2943	pf(k)=pf_i(nlevel+1-k)
2944	print *,'pf=',pf(k)
2945	th(k)=th_i(nlevel+1-k)
2946	t(k)=t_i(nlevel+1-k)
2947	qv(k)=qv_i(nlevel+1-k)
2948	u(k)=u_i(nlevel+1-k)
2949	v(k)=v_i(nlevel+1-k)
2950	hadvt(k,:)=hadvt_i(nlevel+1-k,:)
2951	hadvq(k,:)=hadvq_i(nlevel+1-k,:)
2952	enddo
2953	return
2954	end subroutine read_gabls4
2955	!=====================================================================
2956
2957	! Reads CIRC input files
2958
2959	SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
2960
2961	parameter (ncm_1=49180)
2962	INCLUDE "YOMCST.h"
2963
2964	real albsfc(ncm_1), albsfc_w(ncm_1)
2965	real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
2966	reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
2967	real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
2968	real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
2969	real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
2970	real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
2971	o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
2972	! za= zenital angle
2973	! sza= cosinus angle zenital
2974	real wavn(ncm_1), ssf(ncm_1),za,sza
2975	integer nlev
2976
2977
2978	! Open the files
2979
2980	open (11, file='Tsfc_sza_nlev_case.txt', status='old')
2981	open (12, file='level_input_case.txt', status='old')
2982	open (13, file='layer_input_case.txt', status='old')
2983	open (14, file='aerosol_input_case.txt', status='old')
2984	open (15, file='cloud_input_case.txt', status='old')
2985	open (16, file='sfcalbedo_input_case.txt', status='old')
2986
2987	! Read scalar information
2988	do iskip=1,5
2989	read (11, *)
2990	enddo
2991	read (11, '(i8)') nlev
2992	read (11, '(f10.2)') tsfc
2993	read (11, '(f10.2)') za
2994	read (11, '(f10.4)') sw_dn_toa
2995	sza=cos(za/180.*RPI)
2996	print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
2997	close(11)
2998
2999	! Read level information
3000	read (12, *)
3001	do il=1,nlev
3002	read (12, 302) ilev, z(il), p(il), t(il)
3003	z(il)=z(il)*1000. ! z donne en km
3004	p(il)=p(il)*100. ! p donne en mb
3005	enddo
3006	302 format (i8, f8.3, 2f9.2)
3007	close(12)
3008
3009	! Read layer information (midpoint values)
3010	do iskip=1,3
3011	read (13, *)
3012	enddo
3013	do il=1,nlev-1
3014	read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
3015	n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
3016	f11(il),f12(il)
3017	pm(il)=pm(il)*100.
3018	enddo
3019	303 format (i8, 2f9.2, 10(2x,e13.7))
3020	close(13)
3021
3022	! Read aerosol layer information
3023	do iskip=1,3
3024	read (14, *)
3025	enddo
3026	read (14, '(f10.2)') aer_alpha
3027	read (14, *)
3028	read (14, *)
3029	do il=1,nlev-1
3030	read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
3031	enddo
3032	304 format (i8, f9.5, 2f8.3)
3033	close(14)
3034
3035	! Read cloud information
3036	do iskip=1,3
3037	read (15, *)
3038	enddo
3039	do il=1,nlev-1
3040	read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
3041	lwp(il)=lwp(il)/1000. ! lwp donne en g/kg
3042	iwp(il)=iwp(il)/1000. ! iwp donne en g/kg
3043	reliq(il)=reliq(il)/1000000. ! reliq donne en microns
3044	reice(il)=reice(il)/1000000. ! reice donne en microns
3045	enddo
3046	305 format (i8, f8.3, 4f9.2)
3047	close(15)
3048
3049	! Read surface albedo (weighted & unweighted) and spectral solar irradiance
3050	do iskip=1,6
3051	read (16, *)
3052	enddo
3053	do icm_1=1,ncm_1
3054	read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
3055	enddo
3056	306 format(f10.1, 2f12.5, f14.8)
3057	close(16)
3058
3059	return
3060	end subroutine read_circ
3061	!=====================================================================
3062	! Reads RTMIP input files
3063
3064	SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
3065
3066	INCLUDE "YOMCST.h"
3067
3068	real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
3069	real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
3070	integer nlev
3071
3072
3073	! Open the files
3074
3075	open (11, file='low_resolution_profile.txt', status='old')
3076
3077	! Read level information
3078	read (11, *)
3079	do il=1,nlev_rtmip
3080	read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
3081	enddo
3082	do il=1,nlev_rtmip
3083	play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb
3084	temp(il)=t(nlev_rtmip-il+1)
3085	ovap(il)=h2o(nlev_rtmip-il+1)
3086	oz(il)=o3(nlev_rtmip-il+1)
3087	enddo
3088	do il=1,39
3089	plev(il)=play(il)+(play(il+1)-play(il))/2.
3090	print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
3091	enddo
3092	plev(41)=101300.
3093	302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
3094	close(12)
3095
3096	return
3097	end subroutine read_rtmip
3098	!=====================================================================

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

Download in other formats: