Context Navigation

old_1DUTILS_read_interp.h @ 5447

Last change on this file since 5447 was 5390, checked in by yann meurdesoif, 3 weeks ago

Remove UTF8 character that inihibit fortran parsing with GPU morphosis
Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

File size: 110.0 KB

Line
1	!======================================================================
2	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &
3	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
4	& ,ht_toga,vt_toga,hq_toga,vq_toga)
5	implicit none
6
7	!-------------------------------------------------------------------------
8	! Read TOGA-COARE forcing data
9	!-------------------------------------------------------------------------
10
11	integer nlev_toga,nt_toga
12	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
13	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
14	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
15	real w_toga(nlev_toga,nt_toga)
16	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
17	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
18	character*80 fich_toga
19
20	integer k,ip
21	real bid
22
23	integer iy,im,id,ih
24
25	real plev_min
26
27	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
28
29	open(21,file=trim(fich_toga),form='formatted')
30	read(21,'(a)')
31	do ip = 1, nt_toga
32	read(21,'(a)')
33	read(21,'(a)')
34	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
35	read(21,'(a)')
36	read(21,'(a)')
37
38	do k = 1, nlev_toga
39	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &
40	& ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &
41	& ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
42
43	! conversion in SI units:
44	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
45	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
46	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
47	! no water vapour tendency above 55 hPa
48	if (plev_toga(k,ip) .lt. plev_min) then
49	q_toga(k,ip) = 0.
50	hq_toga(k,ip) = 0.
51	vq_toga(k,ip) =0.
52	endif
53	enddo
54
55	ts_toga(ip)=ts_toga(ip)+273.15 ! K
56	enddo
57	close(21)
58
59	223 format(4i3,6f8.2)
60	230 format(6f9.3,4e11.3)
61
62	return
63	end SUBROUTINE read_togacoare
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 SUBROUTINE read_sandu
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 SUBROUTINE read_astex
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 subroutine catchaxis
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 SUBROUTINE interp_sandu_vertical
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 SUBROUTINE interp_astex_vertical
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 SUBROUTINE read_rico
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 SUBROUTINE interp_sandu_time
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 SUBROUTINE read_armcu
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 SUBROUTINE interp_toga_vertical
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 SUBROUTINE interp_case_vertical
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 SUBROUTINE interp_dice_vertical
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 SUBROUTINE interp_astex_time
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 SUBROUTINE interp_toga_time
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 SUBROUTINE interp_dice_time
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 SUBROUTINE interp_gabls4_time
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 SUBROUTINE interp_armcu_time
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 subroutine readprofiles
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 subroutine readprofile_sandu
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 subroutine readprofile_astex
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 subroutine readprofile_armcu
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 @ 5447

Download in other formats: