Context Navigation

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

Last change on this file since 4586 was 4274, checked in by lguez, 22 months ago

Replace nf_get_var_type by nf90_get_var

The immediate motivation is a bug fix: nf_get_var_type was called
with scalar resul3 or lat, lon, alt, phis instead of array actual
argument for dummy array argument rvals or dvals. Correcting this, we
might as well take the opportunity to use nf90_get_var, so we no
longer need to test NC_DOUBLE and we have half as many calls.

File size: 108.0 KB

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

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

Download in other formats: