Context Navigation

old_1DUTILS_read_interp.h @ 5368

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

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

Context Navigation

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

Download in other formats: