Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

old_1DUTILS_read_interp.h

Last change on this file was 5160, checked in by abarral, 6 months ago
Put .h into modules
File size: 109.3 KB

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

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

Context Navigation

source: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h

Download in other formats: