Context Navigation

source: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h @ 3605

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

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

Download in other formats: