Context Navigation

old_1DUTILS_read_interp.h @ 3679

Last change on this file since 3679 was 3541, checked in by fhourdin, 5 years ago

Gros nettoyage en cours sur le 1D.
Le nouveau lmdz1d.F90 est une coquille vide qui choisit entre
old_lmdz1d.F90 (l'ancien lmdz1d.F90) et scm.F90 (le nouveau au format standard).
Plusieur fichiers sont renommés old_truc, le truc étant au format standard,
nettoyé des anciens formats.
Le 1DUTILS.h est lui même séparé entre des routines génériques venant remplacer
notamment des routines de dyn3d (la vocation d'origine de 1DUTILS.h) et
les routiles de lecture spécifiques mises dans old_1DUTILS.h
On perdra un peu de l'utilister de svn au moment de cette grosse bascule.
Mais les old_ sont faits pour ne plus bouger, et les versions standard
sont en pleine évolution.
Fredho

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.

Context Navigation

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

Download in other formats: