Context Navigation

source: LMDZ5/branches/IPSLCM6.0.10/libf/phylmd/dyn1d/1DUTILS.h @ 5448

Last change on this file since 5448 was 2839, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2785:2838 into testing branch
File size: 173.7 KB

Line
1	#include "conf_gcm.F90"
2
3	!
4	! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
5	!
6	!
7	!
8	SUBROUTINE conf_unicol
9	!
10	#ifdef CPP_IOIPSL
11	use IOIPSL
12	#else
13	! if not using IOIPSL, we still need to use (a local version of) getin
14	use ioipsl_getincom
15	#endif
16	USE print_control_mod, ONLY: lunout
17	IMPLICIT NONE
18	!-----------------------------------------------------------------------
19	! Auteurs : A. Lahellec .
20	!
21	! Declarations :
22	! --------------
23
24	#include "compar1d.h"
25	#include "flux_arp.h"
26	#include "tsoilnudge.h"
27	#include "fcg_gcssold.h"
28	#include "fcg_racmo.h"
29	!
30	!
31	! local:
32	! ------
33
34	! CHARACTER ch172,ch272,ch372,ch412
35
36	!
37	! -------------------------------------------------------------------
38	!
39	! ......... Initilisation parametres du lmdz1D ..........
40	!
41	!---------------------------------------------------------------------
42	! initialisations:
43	! ----------------
44
45	!Config Key = lunout
46	!Config Desc = unite de fichier pour les impressions
47	!Config Def = 6
48	!Config Help = unite de fichier pour les impressions
49	!Config (defaut sortie standard = 6)
50	lunout=6
51	! CALL getin('lunout', lunout)
52	IF (lunout /= 5 .and. lunout /= 6) THEN
53	OPEN(lunout,FILE='lmdz.out')
54	ENDIF
55
56	!Config Key = prt_level
57	!Config Desc = niveau d'impressions de debogage
58	!Config Def = 0
59	!Config Help = Niveau d'impression pour le debogage
60	!Config (0 = minimum d'impression)
61	! prt_level = 0
62	! CALL getin('prt_level',prt_level)
63
64	!-----------------------------------------------------------------------
65	! Parametres de controle du run:
66	!-----------------------------------------------------------------------
67
68	!Config Key = restart
69	!Config Desc = on repart des startphy et start1dyn
70	!Config Def = false
71	!Config Help = les fichiers restart doivent etre renomme en start
72	restart =.false.
73	CALL getin('restart',restart)
74
75	!Config Key = forcing_type
76	!Config Desc = defines the way the SCM is forced:
77	!Config Def = 0
78	!!Config Help = 0 ==> forcing_les = .true.
79	! initial profiles from file prof.inp.001
80	! no forcing by LS convergence ;
81	! surface temperature imposed ;
82	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
83	! = 1 ==> forcing_radconv = .true.
84	! idem forcing_type = 0, but the imposed radiative cooling
85	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
86	! then there is no radiative cooling at all)
87	! = 2 ==> forcing_toga = .true.
88	! initial profiles from TOGA-COARE IFA files
89	! LS convergence and SST imposed from TOGA-COARE IFA files
90	! = 3 ==> forcing_GCM2SCM = .true.
91	! initial profiles from the GCM output
92	! LS convergence imposed from the GCM output
93	! = 4 ==> forcing_twpi = .true.
94	! initial profiles from TWPICE nc files
95	! LS convergence and SST imposed from TWPICE nc files
96	! = 5 ==> forcing_rico = .true.
97	! initial profiles from RICO idealized
98	! LS convergence imposed from RICO (cst)
99	! = 6 ==> forcing_amma = .true.
100	! = 10 ==> forcing_case = .true.
101	! initial profiles from case.nc file
102	! = 40 ==> forcing_GCSSold = .true.
103	! initial profile from GCSS file
104	! LS convergence imposed from GCSS file
105	! = 50 ==> forcing_fire = .true.
106	! = 59 ==> forcing_sandu = .true.
107	! initial profiles from sanduref file: see prof.inp.001
108	! SST varying with time and divergence constante: see ifa_sanduref.txt file
109	! Radiation has to be computed interactively
110	! = 60 ==> forcing_astex = .true.
111	! initial profiles from file: see prof.inp.001
112	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
113	! Radiation has to be computed interactively
114	! = 61 ==> forcing_armcu = .true.
115	! initial profiles from file: see prof.inp.001
116	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
117	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
118	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
119	! Radiation to be switched off
120	! > 100 ==> forcing_case = .true. or forcing_case2 = .true.
121	! initial profiles from case.nc file
122	!
123	forcing_type = 0
124	CALL getin('forcing_type',forcing_type)
125	imp_fcg_gcssold = .false.
126	ts_fcg_gcssold = .false.
127	Tp_fcg_gcssold = .false.
128	Tp_ini_gcssold = .false.
129	xTurb_fcg_gcssold = .false.
130	IF (forcing_type .eq.40) THEN
131	CALL getin('imp_fcg',imp_fcg_gcssold)
132	CALL getin('ts_fcg',ts_fcg_gcssold)
133	CALL getin('tp_fcg',Tp_fcg_gcssold)
134	CALL getin('tp_ini',Tp_ini_gcssold)
135	CALL getin('turb_fcg',xTurb_fcg_gcssold)
136	ENDIF
137
138	!Parametres de forcage
139	!Config Key = tend_t
140	!Config Desc = forcage ou non par advection de T
141	!Config Def = false
142	!Config Help = forcage ou non par advection de T
143	tend_t =0
144	CALL getin('tend_t',tend_t)
145
146	!Config Key = tend_q
147	!Config Desc = forcage ou non par advection de q
148	!Config Def = false
149	!Config Help = forcage ou non par advection de q
150	tend_q =0
151	CALL getin('tend_q',tend_q)
152
153	!Config Key = tend_u
154	!Config Desc = forcage ou non par advection de u
155	!Config Def = false
156	!Config Help = forcage ou non par advection de u
157	tend_u =0
158	CALL getin('tend_u',tend_u)
159
160	!Config Key = tend_v
161	!Config Desc = forcage ou non par advection de v
162	!Config Def = false
163	!Config Help = forcage ou non par advection de v
164	tend_v =0
165	CALL getin('tend_v',tend_v)
166
167	!Config Key = tend_w
168	!Config Desc = forcage ou non par vitesse verticale
169	!Config Def = false
170	!Config Help = forcage ou non par vitesse verticale
171	tend_w =0
172	CALL getin('tend_w',tend_w)
173
174	!Config Key = tend_rayo
175	!Config Desc = forcage ou non par dtrad
176	!Config Def = false
177	!Config Help = forcage ou non par dtrad
178	tend_rayo =0
179	CALL getin('tend_rayo',tend_rayo)
180
181
182	!Config Key = nudge_t
183	!Config Desc = constante de nudging de T
184	!Config Def = false
185	!Config Help = constante de nudging de T
186	nudge_t =0.
187	CALL getin('nudge_t',nudge_t)
188
189	!Config Key = nudge_q
190	!Config Desc = constante de nudging de q
191	!Config Def = false
192	!Config Help = constante de nudging de q
193	nudge_q =0.
194	CALL getin('nudge_q',nudge_q)
195
196	!Config Key = nudge_u
197	!Config Desc = constante de nudging de u
198	!Config Def = false
199	!Config Help = constante de nudging de u
200	nudge_u =0.
201	CALL getin('nudge_u',nudge_u)
202
203	!Config Key = nudge_v
204	!Config Desc = constante de nudging de v
205	!Config Def = false
206	!Config Help = constante de nudging de v
207	nudge_v =0.
208	CALL getin('nudge_v',nudge_v)
209
210	!Config Key = nudge_w
211	!Config Desc = constante de nudging de w
212	!Config Def = false
213	!Config Help = constante de nudging de w
214	nudge_w =0.
215	CALL getin('nudge_w',nudge_w)
216
217
218	!Config Key = iflag_nudge
219	!Config Desc = atmospheric nudging ttype (decimal code)
220	!Config Def = 0
221	!Config Help = 0 ==> no nudging
222	! If digit number n of iflag_nudge is set, then nudging of type n is on
223	! If digit number n of iflag_nudge is not set, then nudging of type n is off
224	! (digits are numbered from the right)
225	iflag_nudge = 0
226	CALL getin('iflag_nudge',iflag_nudge)
227
228	!Config Key = ok_flux_surf
229	!Config Desc = forcage ou non par les flux de surface
230	!Config Def = false
231	!Config Help = forcage ou non par les flux de surface
232	ok_flux_surf =.false.
233	CALL getin('ok_flux_surf',ok_flux_surf)
234
235	!Config Key = ok_prescr_ust
236	!Config Desc = ustar impose ou non
237	!Config Def = false
238	!Config Help = ustar impose ou non
239	ok_prescr_ust = .false.
240	CALL getin('ok_prescr_ust',ok_prescr_ust)
241
242	!Config Key = ok_old_disvert
243	!Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
244	!Config Def = false
245	!Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
246	ok_old_disvert = .FALSE.
247	CALL getin('ok_old_disvert',ok_old_disvert)
248
249	!Config Key = time_ini
250	!Config Desc = meaningless in this case
251	!Config Def = 0.
252	!Config Help =
253	tsurf = 0.
254	CALL getin('time_ini',time_ini)
255
256	!Config Key = rlat et rlon
257	!Config Desc = latitude et longitude
258	!Config Def = 0.0 0.0
259	!Config Help = fixe la position de la colonne
260	xlat = 0.
261	xlon = 0.
262	CALL getin('rlat',xlat)
263	CALL getin('rlon',xlon)
264
265	!Config Key = airephy
266	!Config Desc = Grid cell area
267	!Config Def = 1.e11
268	!Config Help =
269	airefi = 1.e11
270	CALL getin('airephy',airefi)
271
272	!Config Key = nat_surf
273	!Config Desc = surface type
274	!Config Def = 0 (ocean)
275	!Config Help = 0=ocean,1=land,2=glacier,3=banquise
276	nat_surf = 0.
277	CALL getin('nat_surf',nat_surf)
278
279	!Config Key = tsurf
280	!Config Desc = surface temperature
281	!Config Def = 290.
282	!Config Help = not used if type_ts_forcing=1 in lmdz1d.F
283	tsurf = 290.
284	CALL getin('tsurf',tsurf)
285
286	!Config Key = psurf
287	!Config Desc = surface pressure
288	!Config Def = 102400.
289	!Config Help =
290	psurf = 102400.
291	CALL getin('psurf',psurf)
292
293	!Config Key = zsurf
294	!Config Desc = surface altitude
295	!Config Def = 0.
296	!Config Help =
297	zsurf = 0.
298	CALL getin('zsurf',zsurf)
299
300	!Config Key = rugos
301	!Config Desc = coefficient de frottement
302	!Config Def = 0.0001
303	!Config Help = calcul du Cdrag
304	rugos = 0.0001
305	CALL getin('rugos',rugos)
306
307	!Config Key = rugosh
308	!Config Desc = coefficient de frottement
309	!Config Def = rugos
310	!Config Help = calcul du Cdrag
311	rugosh = rugos
312	CALL getin('rugosh',rugosh)
313
314
315
316	!Config Key = snowmass
317	!Config Desc = mass de neige de la surface en kg/m2
318	!Config Def = 0.0000
319	!Config Help = snowmass
320	snowmass = 0.0000
321	CALL getin('snowmass',snowmass)
322
323	!Config Key = wtsurf et wqsurf
324	!Config Desc = ???
325	!Config Def = 0.0 0.0
326	!Config Help =
327	wtsurf = 0.0
328	wqsurf = 0.0
329	CALL getin('wtsurf',wtsurf)
330	CALL getin('wqsurf',wqsurf)
331
332	!Config Key = albedo
333	!Config Desc = albedo
334	!Config Def = 0.09
335	!Config Help =
336	albedo = 0.09
337	CALL getin('albedo',albedo)
338
339	!Config Key = agesno
340	!Config Desc = age de la neige
341	!Config Def = 30.0
342	!Config Help =
343	xagesno = 30.0
344	CALL getin('agesno',xagesno)
345
346	!Config Key = restart_runoff
347	!Config Desc = age de la neige
348	!Config Def = 30.0
349	!Config Help =
350	restart_runoff = 0.0
351	CALL getin('restart_runoff',restart_runoff)
352
353	!Config Key = qsolinp
354	!Config Desc = initial bucket water content (kg/m2) when land (5std)
355	!Config Def = 30.0
356	!Config Help =
357	qsolinp = 1.
358	CALL getin('qsolinp',qsolinp)
359
360	!Config Key = zpicinp
361	!Config Desc = denivellation orographie
362	!Config Def = 0.
363	!Config Help = input brise
364	zpicinp = 0.
365	CALL getin('zpicinp',zpicinp)
366	!Config key = nudge_tsoil
367	!Config Desc = activation of soil temperature nudging
368	!Config Def = .FALSE.
369	!Config Help = ...
370
371	nudge_tsoil=.FALSE.
372	CALL getin('nudge_tsoil',nudge_tsoil)
373
374	!Config key = isoil_nudge
375	!Config Desc = level number where soil temperature is nudged
376	!Config Def = 3
377	!Config Help = ...
378
379	isoil_nudge=3
380	CALL getin('isoil_nudge',isoil_nudge)
381
382	!Config key = Tsoil_nudge
383	!Config Desc = target temperature for tsoil(isoil_nudge)
384	!Config Def = 300.
385	!Config Help = ...
386
387	Tsoil_nudge=300.
388	CALL getin('Tsoil_nudge',Tsoil_nudge)
389
390	!Config key = tau_soil_nudge
391	!Config Desc = nudging relaxation time for tsoil
392	!Config Def = 3600.
393	!Config Help = ...
394
395	tau_soil_nudge=3600.
396	CALL getin('tau_soil_nudge',tau_soil_nudge)
397
398	!----------------------------------------------------------
399	! Param??tres de for??age pour les forcages communs:
400	! Pour les forcages communs: ces entiers valent 0 ou 1
401	! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
402	! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale
403	! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv)
404	! forcages en omega, w, vent geostrophique ou ustar
405	! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging
406	!----------------------------------------------------------
407
408	!Config Key = tadv
409	!Config Desc = forcage ou non par advection totale de T
410	!Config Def = false
411	!Config Help = forcage ou non par advection totale de T
412	tadv =0
413	CALL getin('tadv',tadv)
414
415	!Config Key = tadvv
416	!Config Desc = forcage ou non par advection verticale de T
417	!Config Def = false
418	!Config Help = forcage ou non par advection verticale de T
419	tadvv =0
420	CALL getin('tadvv',tadvv)
421
422	!Config Key = tadvh
423	!Config Desc = forcage ou non par advection horizontale de T
424	!Config Def = false
425	!Config Help = forcage ou non par advection horizontale de T
426	tadvh =0
427	CALL getin('tadvh',tadvh)
428
429	!Config Key = thadv
430	!Config Desc = forcage ou non par advection totale de Theta
431	!Config Def = false
432	!Config Help = forcage ou non par advection totale de Theta
433	thadv =0
434	CALL getin('thadv',thadv)
435
436	!Config Key = thadvv
437	!Config Desc = forcage ou non par advection verticale de Theta
438	!Config Def = false
439	!Config Help = forcage ou non par advection verticale de Theta
440	thadvv =0
441	CALL getin('thadvv',thadvv)
442
443	!Config Key = thadvh
444	!Config Desc = forcage ou non par advection horizontale de Theta
445	!Config Def = false
446	!Config Help = forcage ou non par advection horizontale de Theta
447	thadvh =0
448	CALL getin('thadvh',thadvh)
449
450	!Config Key = qadv
451	!Config Desc = forcage ou non par advection totale de Q
452	!Config Def = false
453	!Config Help = forcage ou non par advection totale de Q
454	qadv =0
455	CALL getin('qadv',qadv)
456
457	!Config Key = qadvv
458	!Config Desc = forcage ou non par advection verticale de Q
459	!Config Def = false
460	!Config Help = forcage ou non par advection verticale de Q
461	qadvv =0
462	CALL getin('qadvv',qadvv)
463
464	!Config Key = qadvh
465	!Config Desc = forcage ou non par advection horizontale de Q
466	!Config Def = false
467	!Config Help = forcage ou non par advection horizontale de Q
468	qadvh =0
469	CALL getin('qadvh',qadvh)
470
471	!Config Key = trad
472	!Config Desc = forcage ou non par tendance radiative
473	!Config Def = false
474	!Config Help = forcage ou non par tendance radiative
475	trad =0
476	CALL getin('trad',trad)
477
478	!Config Key = forc_omega
479	!Config Desc = forcage ou non par omega
480	!Config Def = false
481	!Config Help = forcage ou non par omega
482	forc_omega =0
483	CALL getin('forc_omega',forc_omega)
484
485	!Config Key = forc_w
486	!Config Desc = forcage ou non par w
487	!Config Def = false
488	!Config Help = forcage ou non par w
489	forc_w =0
490	CALL getin('forc_w',forc_w)
491
492	!Config Key = forc_geo
493	!Config Desc = forcage ou non par geo
494	!Config Def = false
495	!Config Help = forcage ou non par geo
496	forc_geo =0
497	CALL getin('forc_geo',forc_geo)
498
499	! Meme chose que ok_precr_ust
500	!Config Key = forc_ustar
501	!Config Desc = forcage ou non par ustar
502	!Config Def = false
503	!Config Help = forcage ou non par ustar
504	forc_ustar =0
505	CALL getin('forc_ustar',forc_ustar)
506	IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
507
508	!Config Key = nudging_u
509	!Config Desc = forcage ou non par nudging sur u
510	!Config Def = false
511	!Config Help = forcage ou non par nudging sur u
512	nudging_u =0
513	CALL getin('nudging_u',nudging_u)
514
515	!Config Key = nudging_v
516	!Config Desc = forcage ou non par nudging sur v
517	!Config Def = false
518	!Config Help = forcage ou non par nudging sur v
519	nudging_v =0
520	CALL getin('nudging_v',nudging_v)
521
522	!Config Key = nudging_w
523	!Config Desc = forcage ou non par nudging sur w
524	!Config Def = false
525	!Config Help = forcage ou non par nudging sur w
526	nudging_w =0
527	CALL getin('nudging_w',nudging_w)
528
529	!Config Key = nudging_q
530	!Config Desc = forcage ou non par nudging sur q
531	!Config Def = false
532	!Config Help = forcage ou non par nudging sur q
533	nudging_q =0
534	CALL getin('nudging_q',nudging_q)
535
536	!Config Key = nudging_t
537	!Config Desc = forcage ou non par nudging sur t
538	!Config Def = false
539	!Config Help = forcage ou non par nudging sur t
540	nudging_t =0
541	CALL getin('nudging_t',nudging_t)
542
543
544
545	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
546	write(lunout,*)' Configuration des parametres du gcm1D: '
547	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
548	write(lunout,*)' restart = ', restart
549	write(lunout,*)' forcing_type = ', forcing_type
550	write(lunout,*)' time_ini = ', time_ini
551	write(lunout,*)' rlat = ', xlat
552	write(lunout,*)' rlon = ', xlon
553	write(lunout,*)' airephy = ', airefi
554	write(lunout,*)' nat_surf = ', nat_surf
555	write(lunout,*)' tsurf = ', tsurf
556	write(lunout,*)' psurf = ', psurf
557	write(lunout,*)' zsurf = ', zsurf
558	write(lunout,*)' rugos = ', rugos
559	write(lunout,*)' snowmass=', snowmass
560	write(lunout,*)' wtsurf = ', wtsurf
561	write(lunout,*)' wqsurf = ', wqsurf
562	write(lunout,*)' albedo = ', albedo
563	write(lunout,*)' xagesno = ', xagesno
564	write(lunout,*)' restart_runoff = ', restart_runoff
565	write(lunout,*)' qsolinp = ', qsolinp
566	write(lunout,*)' zpicinp = ', zpicinp
567	write(lunout,*)' nudge_tsoil = ', nudge_tsoil
568	write(lunout,*)' isoil_nudge = ', isoil_nudge
569	write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
570	write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
571	write(lunout,*)' tadv = ', tadv
572	write(lunout,*)' tadvv = ', tadvv
573	write(lunout,*)' tadvh = ', tadvh
574	write(lunout,*)' thadv = ', thadv
575	write(lunout,*)' thadvv = ', thadvv
576	write(lunout,*)' thadvh = ', thadvh
577	write(lunout,*)' qadv = ', qadv
578	write(lunout,*)' qadvv = ', qadvv
579	write(lunout,*)' qadvh = ', qadvh
580	write(lunout,*)' trad = ', trad
581	write(lunout,*)' forc_omega = ', forc_omega
582	write(lunout,*)' forc_w = ', forc_w
583	write(lunout,*)' forc_geo = ', forc_geo
584	write(lunout,*)' forc_ustar = ', forc_ustar
585	write(lunout,*)' nudging_u = ', nudging_u
586	write(lunout,*)' nudging_v = ', nudging_v
587	write(lunout,*)' nudging_t = ', nudging_t
588	write(lunout,*)' nudging_q = ', nudging_q
589	IF (forcing_type .eq.40) THEN
590	write(lunout,*) '--- Forcing type GCSS Old --- with:'
591	write(lunout,*)'imp_fcg',imp_fcg_gcssold
592	write(lunout,*)'ts_fcg',ts_fcg_gcssold
593	write(lunout,*)'tp_fcg',Tp_fcg_gcssold
594	write(lunout,*)'tp_ini',Tp_ini_gcssold
595	write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
596	ENDIF
597
598	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
599	write(lunout,*)
600	!
601	RETURN
602	END
603	!
604	! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
605	!
606	!
607	SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs, &
608	& ucov,vcov,temp,q,omega2)
609	USE dimphy
610	USE mod_grid_phy_lmdz
611	USE mod_phys_lmdz_para
612	USE iophy
613	USE phys_state_var_mod
614	USE iostart
615	USE write_field_phy
616	USE infotrac
617	use control_mod
618	USE comconst_mod, ONLY: im, jm, lllm
619	USE logic_mod, ONLY: fxyhypb, ysinus
620	USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn
621
622	IMPLICIT NONE
623	!=======================================================
624	! Ecriture du fichier de redemarrage sous format NetCDF
625	!=======================================================
626	! Declarations:
627	! -------------
628	include "dimensions.h"
629	!!#include "control.h"
630	include "netcdf.inc"
631
632	! Arguments:
633	! ----------
634	CHARACTER() fichnom
635	!Al1 plev tronque pour .nc mais plev(klev+1):=0
636	real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev)
637	real :: presnivs(klon,klev)
638	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
639	real :: q(klon,klev,nqtot),omega2(klon,klev)
640	! real :: ug(klev),vg(klev),fcoriolis
641	real :: phis(klon)
642
643	! Variables locales pour NetCDF:
644	! ------------------------------
645	INTEGER iq
646	INTEGER length
647	PARAMETER (length = 100)
648	REAL tab_cntrl(length) ! tableau des parametres du run
649	character*4 nmq(nqtot)
650	character*12 modname
651	character*80 abort_message
652	LOGICAL found
653
654	modname = 'dyn1deta0 : '
655	nmq(1)="vap"
656	nmq(2)="cond"
657	do iq=3,nqtot
658	write(nmq(iq),'("tra",i1)') iq-2
659	enddo
660	print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
661	CALL open_startphy(fichnom)
662	print*,'after open startphy ',fichnom,nmq
663
664	!
665	! Lecture des parametres de controle:
666	!
667	CALL get_var("controle",tab_cntrl)
668
669
670	im = tab_cntrl(1)
671	jm = tab_cntrl(2)
672	lllm = tab_cntrl(3)
673	day_ref = tab_cntrl(4)
674	annee_ref = tab_cntrl(5)
675	! rad = tab_cntrl(6)
676	! omeg = tab_cntrl(7)
677	! g = tab_cntrl(8)
678	! cpp = tab_cntrl(9)
679	! kappa = tab_cntrl(10)
680	! daysec = tab_cntrl(11)
681	! dtvr = tab_cntrl(12)
682	! etot0 = tab_cntrl(13)
683	! ptot0 = tab_cntrl(14)
684	! ztot0 = tab_cntrl(15)
685	! stot0 = tab_cntrl(16)
686	! ang0 = tab_cntrl(17)
687	! pa = tab_cntrl(18)
688	! preff = tab_cntrl(19)
689	!
690	! clon = tab_cntrl(20)
691	! clat = tab_cntrl(21)
692	! grossismx = tab_cntrl(22)
693	! grossismy = tab_cntrl(23)
694	!
695	IF ( tab_cntrl(24).EQ.1. ) THEN
696	fxyhypb =.true.
697	! dzoomx = tab_cntrl(25)
698	! dzoomy = tab_cntrl(26)
699	! taux = tab_cntrl(28)
700	! tauy = tab_cntrl(29)
701	ELSE
702	fxyhypb = .false.
703	ysinus = .false.
704	IF( tab_cntrl(27).EQ.1. ) ysinus =.true.
705	ENDIF
706
707	day_ini = tab_cntrl(30)
708	itau_dyn = tab_cntrl(31)
709	! .................................................................
710	!
711	!
712	! PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
713	!Al1
714	Print*,'day_ref,annee_ref,day_ini,itau_dyn', &
715	& day_ref,annee_ref,day_ini,itau_dyn
716
717	! Lecture des champs
718	!
719	CALL get_field("play",play,found)
720	IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
721	CALL get_field("phi",phi,found)
722	IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
723	CALL get_field("phis",phis,found)
724	IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
725	CALL get_field("presnivs",presnivs,found)
726	IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
727	CALL get_field("ucov",ucov,found)
728	IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
729	CALL get_field("vcov",vcov,found)
730	IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
731	CALL get_field("temp",temp,found)
732	IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
733	CALL get_field("omega2",omega2,found)
734	IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
735	plev(1,klev+1)=0.
736	CALL get_field("plev",plev(:,1:klev),found)
737	IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
738
739	Do iq=1,nqtot
740	CALL get_field("q"//nmq(iq),q(:,:,iq),found)
741	IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
742	EndDo
743
744	CALL close_startphy
745	print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
746	!
747	RETURN
748	END
749	!
750	! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
751	!
752	!
753	SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs, &
754	& ucov,vcov,temp,q,omega2)
755	USE dimphy
756	USE mod_grid_phy_lmdz
757	USE mod_phys_lmdz_para
758	USE phys_state_var_mod
759	USE iostart
760	USE infotrac
761	use control_mod
762	USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
763	USE logic_mod, ONLY: fxyhypb, ysinus
764	USE temps_mod, ONLY: annee_ref,day_end,day_ref,itau_dyn,itaufin
765
766	IMPLICIT NONE
767	!=======================================================
768	! Ecriture du fichier de redemarrage sous format NetCDF
769	!=======================================================
770	! Declarations:
771	! -------------
772	include "dimensions.h"
773	!!#include "control.h"
774	include "netcdf.inc"
775
776	! Arguments:
777	! ----------
778	CHARACTER() fichnom
779	!Al1 plev tronque pour .nc mais plev(klev+1):=0
780	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
781	real :: presnivs(klon,klev)
782	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
783	real :: q(klon,klev,nqtot)
784	real :: omega2(klon,klev),rho(klon,klev+1)
785	! real :: ug(klev),vg(klev),fcoriolis
786	real :: phis(klon)
787
788	! Variables locales pour NetCDF:
789	! ------------------------------
790	INTEGER nid
791	INTEGER ierr
792	INTEGER iq,l
793	INTEGER length
794	PARAMETER (length = 100)
795	REAL tab_cntrl(length) ! tableau des parametres du run
796	character*4 nmq(nqtot)
797	character*20 modname
798	character*80 abort_message
799	!
800	INTEGER nb
801	SAVE nb
802	DATA nb / 0 /
803
804	CALL open_restartphy(fichnom)
805	print*,'redm1 ',fichnom,klon,klev,nqtot
806	nmq(1)="vap"
807	nmq(2)="cond"
808	nmq(3)="tra1"
809	nmq(4)="tra2"
810
811	modname = 'dyn1dredem'
812	ierr = NF_OPEN(fichnom, NF_WRITE, nid)
813	IF (ierr .NE. NF_NOERR) THEN
814	abort_message="Pb. d ouverture "//fichnom
815	CALL abort_gcm('Modele 1D',abort_message,1)
816	ENDIF
817
818	DO l=1,length
819	tab_cntrl(l) = 0.
820	ENDDO
821	tab_cntrl(1) = FLOAT(iim)
822	tab_cntrl(2) = FLOAT(jjm)
823	tab_cntrl(3) = FLOAT(llm)
824	tab_cntrl(4) = FLOAT(day_ref)
825	tab_cntrl(5) = FLOAT(annee_ref)
826	tab_cntrl(6) = rad
827	tab_cntrl(7) = omeg
828	tab_cntrl(8) = g
829	tab_cntrl(9) = cpp
830	tab_cntrl(10) = kappa
831	tab_cntrl(11) = daysec
832	tab_cntrl(12) = dtvr
833	! tab_cntrl(13) = etot0
834	! tab_cntrl(14) = ptot0
835	! tab_cntrl(15) = ztot0
836	! tab_cntrl(16) = stot0
837	! tab_cntrl(17) = ang0
838	! tab_cntrl(18) = pa
839	! tab_cntrl(19) = preff
840	!
841	! ..... parametres pour le zoom ......
842
843	! tab_cntrl(20) = clon
844	! tab_cntrl(21) = clat
845	! tab_cntrl(22) = grossismx
846	! tab_cntrl(23) = grossismy
847	!
848	IF ( fxyhypb ) THEN
849	tab_cntrl(24) = 1.
850	! tab_cntrl(25) = dzoomx
851	! tab_cntrl(26) = dzoomy
852	tab_cntrl(27) = 0.
853	! tab_cntrl(28) = taux
854	! tab_cntrl(29) = tauy
855	ELSE
856	tab_cntrl(24) = 0.
857	! tab_cntrl(25) = dzoomx
858	! tab_cntrl(26) = dzoomy
859	tab_cntrl(27) = 0.
860	tab_cntrl(28) = 0.
861	tab_cntrl(29) = 0.
862	IF( ysinus ) tab_cntrl(27) = 1.
863	ENDIF
864	!Al1 iday_end -> day_end
865	tab_cntrl(30) = FLOAT(day_end)
866	tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
867	!
868	CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
869	!
870
871	! Ecriture/extension de la coordonnee temps
872
873	nb = nb + 1
874
875	! Ecriture des champs
876	!
877	CALL put_field("plev","p interfaces sauf la nulle",plev)
878	CALL put_field("play","",play)
879	CALL put_field("phi","geopotentielle",phi)
880	CALL put_field("phis","geopotentiell de surface",phis)
881	CALL put_field("presnivs","",presnivs)
882	CALL put_field("ucov","",ucov)
883	CALL put_field("vcov","",vcov)
884	CALL put_field("temp","",temp)
885	CALL put_field("omega2","",omega2)
886
887	Do iq=1,nqtot
888	CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs", &
889	& q(:,:,iq))
890	EndDo
891	CALL close_restartphy
892
893	!
894	RETURN
895	END
896	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
897	IMPLICIT NONE
898	!=======================================================================
899	! passage d'un champ de la grille scalaire a la grille physique
900	!=======================================================================
901
902	!-----------------------------------------------------------------------
903	! declarations:
904	! -------------
905
906	INTEGER im,jm,ngrid,nfield
907	REAL pdyn(im,jm,nfield)
908	REAL pfi(ngrid,nfield)
909
910	INTEGER i,j,ifield,ig
911
912	!-----------------------------------------------------------------------
913	! calcul:
914	! -------
915
916	DO ifield=1,nfield
917	! traitement des poles
918	DO i=1,im
919	pdyn(i,1,ifield)=pfi(1,ifield)
920	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
921	ENDDO
922
923	! traitement des point normaux
924	DO j=2,jm-1
925	ig=2+(j-2)*(im-1)
926	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
927	pdyn(im,j,ifield)=pdyn(1,j,ifield)
928	ENDDO
929	ENDDO
930
931	RETURN
932	END
933
934
935
936	SUBROUTINE abort_gcm(modname, message, ierr)
937
938	USE IOIPSL
939	!
940	! Stops the simulation cleanly, closing files and printing various
941	! comments
942	!
943	! Input: modname = name of calling program
944	! message = stuff to print
945	! ierr = severity of situation ( = 0 normal )
946
947	character(len=*) modname
948	integer ierr
949	character(len=*) message
950
951	write(,) 'in abort_gcm'
952	call histclo
953	! call histclo(2)
954	! call histclo(3)
955	! call histclo(4)
956	! call histclo(5)
957	write(,) 'out of histclo'
958	write(,) 'Stopping in ', modname
959	write(,) 'Reason = ',message
960	call getin_dump
961	!
962	if (ierr .eq. 0) then
963	write(,) 'Everything is cool'
964	else
965	write(,) 'Houston, we have a problem ', ierr
966	endif
967	STOP
968	END
969	REAL FUNCTION fq_sat(kelvin, millibar)
970	!
971	IMPLICIT none
972	!======================================================================
973	! Autheur(s): Z.X. Li (LMD/CNRS)
974	! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
975	!======================================================================
976	! Arguments:
977	! kelvin---input-R: temperature en Kelvin
978	! millibar--input-R: pression en mb
979	!
980	! fq_sat----output-R: vapeur d'eau saturante en kg/kg
981	!======================================================================
982	!
983	REAL kelvin, millibar
984	!
985	REAL r2es
986	PARAMETER (r2es=611.14 *18.0153/28.9644)
987	!
988	REAL r3les, r3ies, r3es
989	PARAMETER (R3LES=17.269)
990	PARAMETER (R3IES=21.875)
991	!
992	REAL r4les, r4ies, r4es
993	PARAMETER (R4LES=35.86)
994	PARAMETER (R4IES=7.66)
995	!
996	REAL rtt
997	PARAMETER (rtt=273.16)
998	!
999	REAL retv
1000	PARAMETER (retv=28.9644/18.0153 - 1.0)
1001	!
1002	REAL zqsat
1003	REAL temp, pres
1004	! ------------------------------------------------------------------
1005	!
1006	!
1007	temp = kelvin
1008	pres = millibar * 100.0
1009	! write(,)'kelvin,millibar=',kelvin,millibar
1010	! write(,)'temp,pres=',temp,pres
1011	!
1012	IF (temp .LE. rtt) THEN
1013	r3es = r3ies
1014	r4es = r4ies
1015	ELSE
1016	r3es = r3les
1017	r4es = r4les
1018	ENDIF
1019	!
1020	zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
1021	zqsat=MIN(0.5,ZQSAT)
1022	zqsat=zqsat/(1.-retv *zqsat)
1023	!
1024	fq_sat = zqsat
1025	!
1026	RETURN
1027	END
1028
1029	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
1030	IMPLICIT NONE
1031	!=======================================================================
1032	! passage d'un champ de la grille scalaire a la grille physique
1033	!=======================================================================
1034
1035	!-----------------------------------------------------------------------
1036	! declarations:
1037	! -------------
1038
1039	INTEGER im,jm,ngrid,nfield
1040	REAL pdyn(im,jm,nfield)
1041	REAL pfi(ngrid,nfield)
1042
1043	INTEGER j,ifield,ig
1044
1045	!-----------------------------------------------------------------------
1046	! calcul:
1047	! -------
1048
1049	IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1) &
1050	& STOP 'probleme de dim'
1051	! traitement des poles
1052	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
1053	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
1054
1055	! traitement des point normaux
1056	DO ifield=1,nfield
1057	DO j=2,jm-1
1058	ig=2+(j-2)*(im-1)
1059	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
1060	ENDDO
1061	ENDDO
1062
1063	RETURN
1064	END
1065
1066	SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
1067
1068	! Ancienne version disvert dont on a modifie nom pour utiliser
1069	! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
1070	! (MPL 18092012)
1071	!
1072	! Auteur : P. Le Van .
1073	!
1074	IMPLICIT NONE
1075
1076	include "dimensions.h"
1077	include "paramet.h"
1078	!
1079	!=======================================================================
1080	!
1081	!
1082	! s = sigma ** kappa : coordonnee verticale
1083	! dsig(l) : epaisseur de la couche l ds la coord. s
1084	! sig(l) : sigma a l'interface des couches l et l-1
1085	! ds(l) : distance entre les couches l et l-1 en coord.s
1086	!
1087	!=======================================================================
1088	!
1089	REAL pa,preff
1090	REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
1091	REAL presnivs(llm)
1092	!
1093	! declarations:
1094	! -------------
1095	!
1096	REAL sig(llm+1),dsig(llm)
1097	!
1098	INTEGER l
1099	REAL snorm
1100	REAL alpha,beta,gama,delta,deltaz,h
1101	INTEGER np,ierr
1102	REAL pi,x
1103
1104	!-----------------------------------------------------------------------
1105	!
1106	pi=2.*ASIN(1.)
1107
1108	OPEN(99,file='sigma.def',status='old',form='formatted', &
1109	& iostat=ierr)
1110
1111	!-----------------------------------------------------------------------
1112	! cas 1 on lit les options dans sigma.def:
1113	! ----------------------------------------
1114
1115	IF (ierr.eq.0) THEN
1116
1117	print*,'WARNING!!! on lit les options dans sigma.def'
1118	READ(99,*) deltaz
1119	READ(99,*) h
1120	READ(99,*) beta
1121	READ(99,*) gama
1122	READ(99,*) delta
1123	READ(99,*) np
1124	CLOSE(99)
1125	alpha=deltaz/(llm*h)
1126	!
1127
1128	DO 1 l = 1, llm
1129	dsig(l) = (alpha+(1.-alpha)exp(-beta(llm-l)))* &
1130	& ( (tanh(gamal)/tanh(gamallm))**np + &
1131	& (1.-l/FLOAT(llm))*delta )
1132	1 CONTINUE
1133
1134	sig(1)=1.
1135	DO 101 l=1,llm-1
1136	sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
1137	101 CONTINUE
1138	sig(llm+1)=0.
1139
1140	DO 2 l = 1, llm
1141	dsig(l) = sig(l)-sig(l+1)
1142	2 CONTINUE
1143	!
1144
1145	ELSE
1146	!-----------------------------------------------------------------------
1147	! cas 2 ancienne discretisation (LMD5...):
1148	! ----------------------------------------
1149
1150	PRINT*,'WARNING!!! Ancienne discretisation verticale'
1151
1152	h=7.
1153	snorm = 0.
1154	DO l = 1, llm
1155	x = 2.asin(1.) (FLOAT(l)-0.5) / float(llm+1)
1156	dsig(l) = 1.0 + 7.0 * SIN(x)**2
1157	snorm = snorm + dsig(l)
1158	ENDDO
1159	snorm = 1./snorm
1160	DO l = 1, llm
1161	dsig(l) = dsig(l)*snorm
1162	ENDDO
1163	sig(llm+1) = 0.
1164	DO l = llm, 1, -1
1165	sig(l) = sig(l+1) + dsig(l)
1166	ENDDO
1167
1168	ENDIF
1169
1170
1171	DO l=1,llm
1172	nivsigs(l) = FLOAT(l)
1173	ENDDO
1174
1175	DO l=1,llmp1
1176	nivsig(l)= FLOAT(l)
1177	ENDDO
1178
1179	!
1180	! .... Calculs de ap(l) et de bp(l) ....
1181	! .........................................
1182	!
1183	!
1184	! ..... pa et preff sont lus sur les fichiers start par lectba .....
1185	!
1186
1187	bp(llmp1) = 0.
1188
1189	DO l = 1, llm
1190	!c
1191	!cc ap(l) = 0.
1192	!cc bp(l) = sig(l)
1193
1194	bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
1195	ap(l) = pa * ( sig(l) - bp(l) )
1196	!
1197	ENDDO
1198	ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
1199
1200	PRINT *,' BP '
1201	PRINT *, bp
1202	PRINT *,' AP '
1203	PRINT *, ap
1204
1205	DO l = 1, llm
1206	dpres(l) = bp(l) - bp(l+1)
1207	presnivs(l) = 0.5 ( ap(l)+bp(l)preff + ap(l+1)+bp(l+1)*preff )
1208	ENDDO
1209
1210	PRINT *,' PRESNIVS '
1211	PRINT *,presnivs
1212
1213	RETURN
1214	END
1215
1216	!======================================================================
1217	SUBROUTINE read_tsurf1d(knon,sst_out)
1218
1219	! This subroutine specifies the surface temperature to be used in 1D simulations
1220
1221	USE dimphy, ONLY : klon
1222
1223	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
1224	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
1225
1226	INTEGER :: i
1227	! COMMON defined in lmdz1d.F:
1228	real ts_cur
1229	common /sst_forcing/ts_cur
1230
1231	DO i = 1, knon
1232	sst_out(i) = ts_cur
1233	ENDDO
1234
1235	END SUBROUTINE read_tsurf1d
1236
1237	!===============================================================
1238	subroutine advect_vert(llm,w,dt,q,plev)
1239	!===============================================================
1240	! Schema amont pour l'advection verticale en 1D
1241	! w est la vitesse verticale dp/dt en Pa/s
1242	! Traitement en volumes finis
1243	! d / dt ( zm q ) = delta_z ( omega q )
1244	! d / dt ( zm ) = delta_z ( omega )
1245	! avec zm = delta_z ( p )
1246	! si * designe la valeur au pas de temps t+dt
1247	! zm(l) q(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1248	! zm*(l) -zm(l) = w(l+1) - w(l)
1249	! avec w=omega * dt
1250	!---------------------------------------------------------------
1251	implicit none
1252	! arguments
1253	integer llm
1254	real w(llm+1),q(llm),plev(llm+1),dt
1255
1256	! local
1257	integer l
1258	real zwq(llm+1),zm(llm+1),zw(llm+1)
1259	real qold
1260
1261	!---------------------------------------------------------------
1262
1263	do l=1,llm
1264	zw(l)=dt*w(l)
1265	zm(l)=plev(l)-plev(l+1)
1266	zwq(l)=q(l)*zw(l)
1267	enddo
1268	zwq(llm+1)=0.
1269	zw(llm+1)=0.
1270
1271	do l=1,llm
1272	qold=q(l)
1273	q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1274	print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1275	enddo
1276
1277
1278	return
1279	end
1280
1281	!===============================================================
1282
1283
1284	SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, &
1285	& q,temp,u,v,play)
1286	!itlmd
1287	!----------------------------------------------------------------------
1288	! Calcul de l'advection verticale (ascendance et subsidence) de
1289	! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1290	! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1291	! sans WTG rajouter une advection horizontale
1292	!----------------------------------------------------------------------
1293	implicit none
1294	#include "YOMCST.h"
1295	! argument
1296	integer llm
1297	real omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1298	real d_u_va(llm), d_v_va(llm)
1299	real q(llm,3),temp(llm)
1300	real u(llm),v(llm)
1301	real play(llm)
1302	! interne
1303	integer l
1304	real alpha,omgdown,omgup
1305
1306	do l= 1,llm
1307	if(l.eq.1) then
1308	!si omgup pour la couche 1, alors tendance nulle
1309	omgdown=max(omega(2),0.0)
1310	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1311	d_t_va(l)= alpha(omgdown)-omgdown(temp(l)-temp(l+1)) &
1312	& /(play(l)-play(l+1))
1313
1314	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))
1315
1316	d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))
1317	d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))
1318
1319
1320	elseif(l.eq.llm) then
1321	omgup=min(omega(l),0.0)
1322	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1323	d_t_va(l)= alpha*(omgup)- &
1324
1325	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1326	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1327	d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1328	d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1329	d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1330
1331	else
1332	omgup=min(omega(l),0.0)
1333	omgdown=max(omega(l+1),0.0)
1334	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1335	d_t_va(l)= alpha(omgup+omgdown)-omgdown(temp(l)-temp(l+1)) &
1336	& /(play(l)-play(l+1))- &
1337	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1338	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1339	! print*, ' ??? '
1340
1341	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:)) &
1342	& /(play(l)-play(l+1))- &
1343	& omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1344	d_u_va(l)= -omgdown*(u(l)-u(l+1)) &
1345	& /(play(l)-play(l+1))- &
1346	& omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1347	d_v_va(l)= -omgdown*(v(l)-v(l+1)) &
1348	& /(play(l)-play(l+1))- &
1349	& omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1350
1351	endif
1352
1353	enddo
1354	!fin itlmd
1355	return
1356	end
1357	! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1358	SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va, &
1359	& q,temp,u,v,play)
1360	!itlmd
1361	!----------------------------------------------------------------------
1362	! Calcul de l'advection verticale (ascendance et subsidence) de
1363	! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1364	! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1365	! sans WTG rajouter une advection horizontale
1366	!----------------------------------------------------------------------
1367	implicit none
1368	#include "YOMCST.h"
1369	! argument
1370	integer llm,nqtot
1371	real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1372	! real d_u_va(llm), d_v_va(llm)
1373	real q(llm,nqtot),temp(llm)
1374	real u(llm),v(llm)
1375	real play(llm)
1376	real cor(llm)
1377	! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1378	real dph(llm),dqdp(llm),dtdp(llm)
1379	! interne
1380	integer k
1381	real omdn,omup
1382
1383	! dudp=0.
1384	! dvdp=0.
1385	dqdp=0.
1386	dtdp=0.
1387	! d_u_va=0.
1388	! d_v_va=0.
1389
1390	cor(:) = rkappatemp(1.+q(:,1)rv/rd)/(play(1.+q(:,1)))
1391
1392
1393	do k=2,llm-1
1394
1395	dph (k-1) = (play(k )- play(k-1 ))
1396	! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
1397	! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1)
1398	dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1)
1399	dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1)
1400
1401	enddo
1402
1403	! dudp ( llm ) = dudp ( llm-1 )
1404	! dvdp ( llm ) = dvdp ( llm-1 )
1405	dqdp ( llm ) = dqdp ( llm-1 )
1406	dtdp ( llm ) = dtdp ( llm-1 )
1407
1408	do k=2,llm-1
1409	omdn=max(0.0,omega(k+1))
1410	omup=min(0.0,omega( k ))
1411
1412	! d_u_va(k) = -omdndudp(k)-omupdudp(k-1)
1413	! d_v_va(k) = -omdndvdp(k)-omupdvdp(k-1)
1414	d_q_va(k,1)= -omdndqdp(k)-omupdqdp(k-1)
1415	d_t_va(k) = -omdndtdp(k)-omupdtdp(k-1)+(omup+omdn)*cor(k)
1416	enddo
1417
1418	omdn=max(0.0,omega( 2 ))
1419	omup=min(0.0,omega(llm))
1420	! d_u_va( 1 ) = -omdn*dudp( 1 )
1421	! d_u_va(llm) = -omup*dudp(llm)
1422	! d_v_va( 1 ) = -omdn*dvdp( 1 )
1423	! d_v_va(llm) = -omup*dvdp(llm)
1424	d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1425	d_q_va(llm,1) = -omup*dqdp(llm)
1426	d_t_va( 1 ) = -omdndtdp( 1 )+omdncor( 1 )
1427	d_t_va(llm) = -omupdtdp(llm)!+omupcor(llm)
1428
1429	! if(abs(rlat(1))>10.) then
1430	! Calculate the tendency due agestrophic motions
1431	! du_age = fcoriolis*(v-vg)
1432	! dv_age = fcoriolis*(ug-u)
1433	! endif
1434
1435	! call writefield_phy('d_t_va',d_t_va,llm)
1436
1437	return
1438	end
1439
1440	!======================================================================
1441	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &
1442	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
1443	& ,ht_toga,vt_toga,hq_toga,vq_toga)
1444	implicit none
1445
1446	!-------------------------------------------------------------------------
1447	! Read TOGA-COARE forcing data
1448	!-------------------------------------------------------------------------
1449
1450	integer nlev_toga,nt_toga
1451	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1452	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1453	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1454	real w_toga(nlev_toga,nt_toga)
1455	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1456	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1457	character*80 fich_toga
1458
1459	integer k,ip
1460	real bid
1461
1462	integer iy,im,id,ih
1463
1464	real plev_min
1465
1466	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1467
1468	open(21,file=trim(fich_toga),form='formatted')
1469	read(21,'(a)')
1470	do ip = 1, nt_toga
1471	read(21,'(a)')
1472	read(21,'(a)')
1473	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1474	read(21,'(a)')
1475	read(21,'(a)')
1476
1477	do k = 1, nlev_toga
1478	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &
1479	& ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &
1480	& ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1481
1482	! conversion in SI units:
1483	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
1484	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
1485	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1486	! no water vapour tendency above 55 hPa
1487	if (plev_toga(k,ip) .lt. plev_min) then
1488	q_toga(k,ip) = 0.
1489	hq_toga(k,ip) = 0.
1490	vq_toga(k,ip) =0.
1491	endif
1492	enddo
1493
1494	ts_toga(ip)=ts_toga(ip)+273.15 ! K
1495	enddo
1496	close(21)
1497
1498	223 format(4i3,6f8.2)
1499	230 format(6f9.3,4e11.3)
1500
1501	return
1502	end
1503
1504	!-------------------------------------------------------------------------
1505	SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
1506	implicit none
1507
1508	!-------------------------------------------------------------------------
1509	! Read I.SANDU case forcing data
1510	!-------------------------------------------------------------------------
1511
1512	integer nlev_sandu,nt_sandu
1513	real ts_sandu(nt_sandu)
1514	character*80 fich_sandu
1515
1516	integer ip
1517	integer iy,im,id,ih
1518
1519	real plev_min
1520
1521	print*,'nlev_sandu',nlev_sandu
1522	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1523
1524	open(21,file=trim(fich_sandu),form='formatted')
1525	read(21,'(a)')
1526	do ip = 1, nt_sandu
1527	read(21,'(a)')
1528	read(21,'(a)')
1529	read(21,223) iy, im, id, ih, ts_sandu(ip)
1530	print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
1531	enddo
1532	close(21)
1533
1534	223 format(4i3,f8.2)
1535
1536	return
1537	end
1538
1539	!=====================================================================
1540	!-------------------------------------------------------------------------
1541	SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, &
1542	& ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
1543	implicit none
1544
1545	!-------------------------------------------------------------------------
1546	! Read Astex case forcing data
1547	!-------------------------------------------------------------------------
1548
1549	integer nlev_astex,nt_astex
1550	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
1551	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
1552	character*80 fich_astex
1553
1554	integer ip
1555	integer iy,im,id,ih
1556
1557	real plev_min
1558
1559	print*,'nlev_astex',nlev_astex
1560	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1561
1562	open(21,file=trim(fich_astex),form='formatted')
1563	read(21,'(a)')
1564	read(21,'(a)')
1565	do ip = 1, nt_astex
1566	read(21,'(a)')
1567	read(21,'(a)')
1568	read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), &
1569	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
1570	ts_astex(ip)=ts_astex(ip)+273.15
1571	print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), &
1572	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
1573	enddo
1574	close(21)
1575
1576	223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
1577
1578	return
1579	end
1580	!=====================================================================
1581	subroutine read_twpice(fich_twpice,nlevel,ntime &
1582	& ,T_srf,plev,T,q,u,v,omega &
1583	& ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1584
1585	!program reading forcings of the TWP-ICE experiment
1586
1587	! use netcdf
1588
1589	implicit none
1590
1591	#include "netcdf.inc"
1592
1593	integer ntime,nlevel
1594	integer l,k
1595	character*80 :: fich_twpice
1596	real*8 time(ntime)
1597	real*8 lat, lon, alt, phis
1598	real*8 lev(nlevel)
1599	real*8 plev(nlevel,ntime)
1600
1601	real*8 T(nlevel,ntime)
1602	real*8 q(nlevel,ntime),u(nlevel,ntime)
1603	real*8 v(nlevel,ntime)
1604	real*8 omega(nlevel,ntime), div(nlevel,ntime)
1605	real*8 T_adv_h(nlevel,ntime)
1606	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1607	real*8 q_adv_v(nlevel,ntime)
1608	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1609	real*8 s_adv_v(nlevel,ntime)
1610	real*8 p_srf_aver(ntime), p_srf_center(ntime)
1611	real*8 T_srf(ntime)
1612
1613	integer nid, ierr
1614	integer nbvar3d
1615	parameter(nbvar3d=20)
1616	integer var3didin(nbvar3d)
1617
1618	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1619	if (ierr.NE.NF_NOERR) then
1620	write(,) 'ERROR: Pb opening forcings cdf file '
1621	write(,) NF_STRERROR(ierr)
1622	stop ""
1623	endif
1624
1625	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1626	if(ierr/=NF_NOERR) then
1627	write(,) NF_STRERROR(ierr)
1628	stop 'lat'
1629	endif
1630
1631	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1632	if(ierr/=NF_NOERR) then
1633	write(,) NF_STRERROR(ierr)
1634	stop 'lon'
1635	endif
1636
1637	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1638	if(ierr/=NF_NOERR) then
1639	write(,) NF_STRERROR(ierr)
1640	stop 'alt'
1641	endif
1642
1643	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1644	if(ierr/=NF_NOERR) then
1645	write(,) NF_STRERROR(ierr)
1646	stop 'phis'
1647	endif
1648
1649	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1650	if(ierr/=NF_NOERR) then
1651	write(,) NF_STRERROR(ierr)
1652	stop 'T'
1653	endif
1654
1655	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1656	if(ierr/=NF_NOERR) then
1657	write(,) NF_STRERROR(ierr)
1658	stop 'q'
1659	endif
1660
1661	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1662	if(ierr/=NF_NOERR) then
1663	write(,) NF_STRERROR(ierr)
1664	stop 'u'
1665	endif
1666
1667	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1668	if(ierr/=NF_NOERR) then
1669	write(,) NF_STRERROR(ierr)
1670	stop 'v'
1671	endif
1672
1673	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1674	if(ierr/=NF_NOERR) then
1675	write(,) NF_STRERROR(ierr)
1676	stop 'omega'
1677	endif
1678
1679	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1680	if(ierr/=NF_NOERR) then
1681	write(,) NF_STRERROR(ierr)
1682	stop 'div'
1683	endif
1684
1685	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1686	if(ierr/=NF_NOERR) then
1687	write(,) NF_STRERROR(ierr)
1688	stop 'T_adv_h'
1689	endif
1690
1691	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1692	if(ierr/=NF_NOERR) then
1693	write(,) NF_STRERROR(ierr)
1694	stop 'T_adv_v'
1695	endif
1696
1697	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1698	if(ierr/=NF_NOERR) then
1699	write(,) NF_STRERROR(ierr)
1700	stop 'q_adv_h'
1701	endif
1702
1703	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1704	if(ierr/=NF_NOERR) then
1705	write(,) NF_STRERROR(ierr)
1706	stop 'q_adv_v'
1707	endif
1708
1709	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1710	if(ierr/=NF_NOERR) then
1711	write(,) NF_STRERROR(ierr)
1712	stop 's'
1713	endif
1714
1715	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1716	if(ierr/=NF_NOERR) then
1717	write(,) NF_STRERROR(ierr)
1718	stop 's_adv_h'
1719	endif
1720
1721	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1722	if(ierr/=NF_NOERR) then
1723	write(,) NF_STRERROR(ierr)
1724	stop 's_adv_v'
1725	endif
1726
1727	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1728	if(ierr/=NF_NOERR) then
1729	write(,) NF_STRERROR(ierr)
1730	stop 'p_srf_aver'
1731	endif
1732
1733	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1734	if(ierr/=NF_NOERR) then
1735	write(,) NF_STRERROR(ierr)
1736	stop 'p_srf_center'
1737	endif
1738
1739	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1740	if(ierr/=NF_NOERR) then
1741	write(,) NF_STRERROR(ierr)
1742	stop 'T_srf'
1743	endif
1744
1745	!dimensions lecture
1746	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1747
1748	!pressure
1749	do l=1,ntime
1750	do k=1,nlevel
1751	plev(k,l)=lev(k)
1752	enddo
1753	enddo
1754
1755	#ifdef NC_DOUBLE
1756	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1757	#else
1758	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1759	#endif
1760	if(ierr/=NF_NOERR) then
1761	write(,) NF_STRERROR(ierr)
1762	stop "getvarup"
1763	endif
1764	! write(,)'lecture lat ok',lat
1765
1766	#ifdef NC_DOUBLE
1767	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1768	#else
1769	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1770	#endif
1771	if(ierr/=NF_NOERR) then
1772	write(,) NF_STRERROR(ierr)
1773	stop "getvarup"
1774	endif
1775	! write(,)'lecture lon ok',lon
1776
1777	#ifdef NC_DOUBLE
1778	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1779	#else
1780	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1781	#endif
1782	if(ierr/=NF_NOERR) then
1783	write(,) NF_STRERROR(ierr)
1784	stop "getvarup"
1785	endif
1786	! write(,)'lecture alt ok',alt
1787
1788	#ifdef NC_DOUBLE
1789	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1790	#else
1791	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1792	#endif
1793	if(ierr/=NF_NOERR) then
1794	write(,) NF_STRERROR(ierr)
1795	stop "getvarup"
1796	endif
1797	! write(,)'lecture phis ok',phis
1798
1799	#ifdef NC_DOUBLE
1800	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1801	#else
1802	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1803	#endif
1804	if(ierr/=NF_NOERR) then
1805	write(,) NF_STRERROR(ierr)
1806	stop "getvarup"
1807	endif
1808	! write(,)'lecture T ok'
1809
1810	#ifdef NC_DOUBLE
1811	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1812	#else
1813	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1814	#endif
1815	if(ierr/=NF_NOERR) then
1816	write(,) NF_STRERROR(ierr)
1817	stop "getvarup"
1818	endif
1819	! write(,)'lecture q ok'
1820	!q in kg/kg
1821	do l=1,ntime
1822	do k=1,nlevel
1823	q(k,l)=q(k,l)/1000.
1824	enddo
1825	enddo
1826	#ifdef NC_DOUBLE
1827	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1828	#else
1829	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1830	#endif
1831	if(ierr/=NF_NOERR) then
1832	write(,) NF_STRERROR(ierr)
1833	stop "getvarup"
1834	endif
1835	! write(,)'lecture u ok'
1836
1837	#ifdef NC_DOUBLE
1838	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1839	#else
1840	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1841	#endif
1842	if(ierr/=NF_NOERR) then
1843	write(,) NF_STRERROR(ierr)
1844	stop "getvarup"
1845	endif
1846	! write(,)'lecture v ok'
1847
1848	#ifdef NC_DOUBLE
1849	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1850	#else
1851	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1852	#endif
1853	if(ierr/=NF_NOERR) then
1854	write(,) NF_STRERROR(ierr)
1855	stop "getvarup"
1856	endif
1857	! write(,)'lecture omega ok'
1858	!omega in mb/hour
1859	do l=1,ntime
1860	do k=1,nlevel
1861	omega(k,l)=omega(k,l)*100./3600.
1862	enddo
1863	enddo
1864
1865	#ifdef NC_DOUBLE
1866	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1867	#else
1868	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1869	#endif
1870	if(ierr/=NF_NOERR) then
1871	write(,) NF_STRERROR(ierr)
1872	stop "getvarup"
1873	endif
1874	! write(,)'lecture div ok'
1875
1876	#ifdef NC_DOUBLE
1877	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1878	#else
1879	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1880	#endif
1881	if(ierr/=NF_NOERR) then
1882	write(,) NF_STRERROR(ierr)
1883	stop "getvarup"
1884	endif
1885	! write(,)'lecture T_adv_h ok'
1886	!T adv in K/s
1887	do l=1,ntime
1888	do k=1,nlevel
1889	T_adv_h(k,l)=T_adv_h(k,l)/3600.
1890	enddo
1891	enddo
1892
1893
1894	#ifdef NC_DOUBLE
1895	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1896	#else
1897	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1898	#endif
1899	if(ierr/=NF_NOERR) then
1900	write(,) NF_STRERROR(ierr)
1901	stop "getvarup"
1902	endif
1903	! write(,)'lecture T_adv_v ok'
1904	!T adv in K/s
1905	do l=1,ntime
1906	do k=1,nlevel
1907	T_adv_v(k,l)=T_adv_v(k,l)/3600.
1908	enddo
1909	enddo
1910
1911	#ifdef NC_DOUBLE
1912	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1913	#else
1914	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1915	#endif
1916	if(ierr/=NF_NOERR) then
1917	write(,) NF_STRERROR(ierr)
1918	stop "getvarup"
1919	endif
1920	! write(,)'lecture q_adv_h ok'
1921	!q adv in kg/kg/s
1922	do l=1,ntime
1923	do k=1,nlevel
1924	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1925	enddo
1926	enddo
1927
1928
1929	#ifdef NC_DOUBLE
1930	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1931	#else
1932	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1933	#endif
1934	if(ierr/=NF_NOERR) then
1935	write(,) NF_STRERROR(ierr)
1936	stop "getvarup"
1937	endif
1938	! write(,)'lecture q_adv_v ok'
1939	!q adv in kg/kg/s
1940	do l=1,ntime
1941	do k=1,nlevel
1942	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1943	enddo
1944	enddo
1945
1946
1947	#ifdef NC_DOUBLE
1948	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1949	#else
1950	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1951	#endif
1952	if(ierr/=NF_NOERR) then
1953	write(,) NF_STRERROR(ierr)
1954	stop "getvarup"
1955	endif
1956
1957	#ifdef NC_DOUBLE
1958	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1959	#else
1960	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1961	#endif
1962	if(ierr/=NF_NOERR) then
1963	write(,) NF_STRERROR(ierr)
1964	stop "getvarup"
1965	endif
1966
1967	#ifdef NC_DOUBLE
1968	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1969	#else
1970	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1971	#endif
1972	if(ierr/=NF_NOERR) then
1973	write(,) NF_STRERROR(ierr)
1974	stop "getvarup"
1975	endif
1976
1977	#ifdef NC_DOUBLE
1978	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1979	#else
1980	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1981	#endif
1982	if(ierr/=NF_NOERR) then
1983	write(,) NF_STRERROR(ierr)
1984	stop "getvarup"
1985	endif
1986
1987	#ifdef NC_DOUBLE
1988	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1989	#else
1990	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1991	#endif
1992	if(ierr/=NF_NOERR) then
1993	write(,) NF_STRERROR(ierr)
1994	stop "getvarup"
1995	endif
1996
1997	#ifdef NC_DOUBLE
1998	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1999	#else
2000	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
2001	#endif
2002	if(ierr/=NF_NOERR) then
2003	write(,) NF_STRERROR(ierr)
2004	stop "getvarup"
2005	endif
2006	! write(,)'lecture T_srf ok', T_srf
2007
2008	return
2009	end subroutine read_twpice
2010	!=====================================================================
2011	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
2012
2013	! use netcdf
2014
2015	implicit none
2016	#include "netcdf.inc"
2017	integer nid,ttm,llm
2018	real*8 time(ttm)
2019	real*8 lev(llm)
2020	integer ierr
2021
2022	integer timevar,levvar
2023	integer timelen,levlen
2024	integer timedimin,levdimin
2025
2026	! Control & lecture on dimensions
2027	! ===============================
2028	ierr=NF_INQ_DIMID(nid,"time",timedimin)
2029	ierr=NF_INQ_VARID(nid,"time",timevar)
2030	if (ierr.NE.NF_NOERR) then
2031	write(,) 'ERROR: Field <time> is missing'
2032	stop ""
2033	endif
2034	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
2035
2036	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
2037	ierr=NF_INQ_VARID(nid,"lev",levvar)
2038	if (ierr.NE.NF_NOERR) then
2039	write(,) 'ERROR: Field <lev> is lacking'
2040	stop ""
2041	endif
2042	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
2043
2044	if((timelen/=ttm).or.(levlen/=llm)) then
2045	write(,) 'ERROR: Not the good lenght for axis'
2046	write(,) 'longitude: ',timelen,ttm+1
2047	write(,) 'latitude: ',levlen,llm
2048	stop ""
2049	endif
2050
2051	!#ifdef NC_DOUBLE
2052	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
2053	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
2054	!#else
2055	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
2056	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
2057	!#endif
2058
2059	return
2060	end
2061	!=====================================================================
2062
2063	SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof &
2064	& ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof &
2065	& ,omega_prof,o3mmr_prof &
2066	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
2067	& ,omega_mod,o3mmr_mod,mxcalc)
2068
2069	implicit none
2070
2071	#include "dimensions.h"
2072
2073	!-------------------------------------------------------------------------
2074	! Vertical interpolation of SANDUREF forcing data onto model levels
2075	!-------------------------------------------------------------------------
2076
2077	integer nlevmax
2078	parameter (nlevmax=41)
2079	integer nlev_sandu,mxcalc
2080	! real play(llm), plev_prof(nlevmax)
2081	! real t_prof(nlevmax),q_prof(nlevmax)
2082	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2083	! real ht_prof(nlevmax),vt_prof(nlevmax)
2084	! real hq_prof(nlevmax),vq_prof(nlevmax)
2085
2086	real play(llm), plev_prof(nlev_sandu)
2087	real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
2088	real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
2089	real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
2090
2091	real t_mod(llm),thl_mod(llm),q_mod(llm)
2092	real u_mod(llm),v_mod(llm), w_mod(llm)
2093	real omega_mod(llm),o3mmr_mod(llm)
2094
2095	integer l,k,k1,k2
2096	real frac,frac1,frac2,fact
2097
2098	do l = 1, llm
2099
2100	if (play(l).ge.plev_prof(nlev_sandu)) then
2101
2102	mxcalc=l
2103	k1=0
2104	k2=0
2105
2106	if (play(l).le.plev_prof(1)) then
2107
2108	do k = 1, nlev_sandu-1
2109	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2110	k1=k
2111	k2=k+1
2112	endif
2113	enddo
2114
2115	if (k1.eq.0 .or. k2.eq.0) then
2116	write(,) 'PB! k1, k2 = ',k1,k2
2117	write(,) 'l,play(l) = ',l,play(l)/100
2118	do k = 1, nlev_sandu-1
2119	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2120	enddo
2121	endif
2122
2123	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2124	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2125	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
2126	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2127	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2128	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2129	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2130	omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
2131	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
2132
2133	else !play>plev_prof(1)
2134
2135	k1=1
2136	k2=2
2137	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2138	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2139	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2140	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
2141	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
2142	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2143	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2144	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2145	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
2146	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
2147
2148	endif ! play.le.plev_prof(1)
2149
2150	else ! above max altitude of forcing file
2151
2152	!jyg
2153	fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
2154	fact = max(fact,0.) !jyg
2155	fact = exp(-fact) !jyg
2156	t_mod(l)= t_prof(nlev_sandu) !jyg
2157	thl_mod(l)= thl_prof(nlev_sandu) !jyg
2158	q_mod(l)= q_prof(nlev_sandu)*fact !jyg
2159	u_mod(l)= u_prof(nlev_sandu)*fact !jyg
2160	v_mod(l)= v_prof(nlev_sandu)*fact !jyg
2161	w_mod(l)= w_prof(nlev_sandu)*fact !jyg
2162	omega_mod(l)= omega_prof(nlev_sandu)*fact !jyg
2163	o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact !jyg
2164
2165	endif ! play
2166
2167	enddo ! l
2168
2169	do l = 1,llm
2170	! print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
2171	! $ l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
2172	enddo
2173
2174	return
2175	end
2176	!=====================================================================
2177	SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof &
2178	& ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof &
2179	& ,w_prof,tke_prof,o3mmr_prof &
2180	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
2181	& ,tke_mod,o3mmr_mod,mxcalc)
2182
2183	implicit none
2184
2185	#include "dimensions.h"
2186
2187	!-------------------------------------------------------------------------
2188	! Vertical interpolation of Astex forcing data onto model levels
2189	!-------------------------------------------------------------------------
2190
2191	integer nlevmax
2192	parameter (nlevmax=41)
2193	integer nlev_astex,mxcalc
2194	! real play(llm), plev_prof(nlevmax)
2195	! real t_prof(nlevmax),qv_prof(nlevmax)
2196	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2197	! real ht_prof(nlevmax),vt_prof(nlevmax)
2198	! real hq_prof(nlevmax),vq_prof(nlevmax)
2199
2200	real play(llm), plev_prof(nlev_astex)
2201	real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
2202	real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
2203	real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
2204	real qt_prof(nlev_astex),tke_prof(nlev_astex)
2205
2206	real t_mod(llm),thl_mod(llm),qv_mod(llm)
2207	real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
2208	real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
2209
2210	integer l,k,k1,k2
2211	real frac,frac1,frac2,fact
2212
2213	do l = 1, llm
2214
2215	if (play(l).ge.plev_prof(nlev_astex)) then
2216
2217	mxcalc=l
2218	k1=0
2219	k2=0
2220
2221	if (play(l).le.plev_prof(1)) then
2222
2223	do k = 1, nlev_astex-1
2224	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2225	k1=k
2226	k2=k+1
2227	endif
2228	enddo
2229
2230	if (k1.eq.0 .or. k2.eq.0) then
2231	write(,) 'PB! k1, k2 = ',k1,k2
2232	write(,) 'l,play(l) = ',l,play(l)/100
2233	do k = 1, nlev_astex-1
2234	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2235	enddo
2236	endif
2237
2238	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2239	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2240	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
2241	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
2242	ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
2243	qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
2244	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2245	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2246	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2247	tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
2248	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
2249
2250	else !play>plev_prof(1)
2251
2252	k1=1
2253	k2=2
2254	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2255	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2256	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2257	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
2258	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
2259	ql_mod(l)= frac1ql_prof(k1) - frac2ql_prof(k2)
2260	qt_mod(l)= frac1qt_prof(k1) - frac2qt_prof(k2)
2261	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2262	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2263	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2264	tke_mod(l)= frac1tke_prof(k1) - frac2tke_prof(k2)
2265	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
2266
2267	endif ! play.le.plev_prof(1)
2268
2269	else ! above max altitude of forcing file
2270
2271	!jyg
2272	fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
2273	fact = max(fact,0.) !jyg
2274	fact = exp(-fact) !jyg
2275	t_mod(l)= t_prof(nlev_astex) !jyg
2276	thl_mod(l)= thl_prof(nlev_astex) !jyg
2277	qv_mod(l)= qv_prof(nlev_astex)*fact !jyg
2278	ql_mod(l)= ql_prof(nlev_astex)*fact !jyg
2279	qt_mod(l)= qt_prof(nlev_astex)*fact !jyg
2280	u_mod(l)= u_prof(nlev_astex)*fact !jyg
2281	v_mod(l)= v_prof(nlev_astex)*fact !jyg
2282	w_mod(l)= w_prof(nlev_astex)*fact !jyg
2283	tke_mod(l)= tke_prof(nlev_astex)*fact !jyg
2284	o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact !jyg
2285
2286	endif ! play
2287
2288	enddo ! l
2289
2290	do l = 1,llm
2291	! print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
2292	! $ l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
2293	enddo
2294
2295	return
2296	end
2297
2298	!======================================================================
2299	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play &
2300	& ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &
2301	& ,dth_dyn,dqh_dyn)
2302	implicit none
2303
2304	!-------------------------------------------------------------------------
2305	! Read RICO forcing data
2306	!-------------------------------------------------------------------------
2307	#include "dimensions.h"
2308
2309
2310	integer nlev_rico
2311	real ts_rico,ps_rico
2312	real t_rico(llm),q_rico(llm)
2313	real u_rico(llm),v_rico(llm)
2314	real w_rico(llm)
2315	real dth_dyn(llm)
2316	real dqh_dyn(llm)
2317
2318
2319	real play(llm),zlay(llm)
2320
2321
2322	real prico(nlev_rico),zrico(nlev_rico)
2323
2324	character*80 fich_rico
2325
2326	integer k,l
2327
2328
2329	print*,fich_rico
2330	open(21,file=trim(fich_rico),form='formatted')
2331	do k=1,llm
2332	zlay(k)=0.
2333	enddo
2334
2335	read(21,*) ps_rico,ts_rico
2336	prico(1)=ps_rico
2337	zrico(1)=0.0
2338	do l=2,nlev_rico
2339	read(21,*) k,prico(l),zrico(l)
2340	enddo
2341	close(21)
2342
2343	do k=1,llm
2344	do l=1,80
2345	if(prico(l)>play(k)) then
2346	if(play(k)>prico(l+1)) then
2347	zlay(k)=zrico(l)+(play(k)-prico(l)) * &
2348	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2349	else
2350	zlay(k)=zrico(l)+(play(k)-prico(80))* &
2351	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
2352	endif
2353	endif
2354	enddo
2355	print*,k,zlay(k)
2356	! U
2357	if(0 < zlay(k) .and. zlay(k) < 4000) then
2358	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2359	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2360	u_rico(k)= -1.9 + (30.0 + 1.9) / &
2361	& (12000 - 4000) * (zlay(k) - 4000)
2362	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2363	u_rico(k)=30.0
2364	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2365	u_rico(k)=30.0 - (30.0) / &
2366	& (20000 - 13000) * (zlay(k) - 13000)
2367	else
2368	u_rico(k)=0.0
2369	endif
2370
2371	!Q_v
2372	if(0 < zlay(k) .and. zlay(k) < 740) then
2373	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2374	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2375	q_rico(k)=13.8 + (2.4 - 13.8) / &
2376	& (3260 - 740) * (zlay(k) - 740)
2377	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2378	q_rico(k)=2.4 + (1.8 - 2.4) / &
2379	& (4000 - 3260) * (zlay(k) - 3260)
2380	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2381	q_rico(k)=1.8 + (0 - 1.8) / &
2382	& (9000 - 4000) * (zlay(k) - 4000)
2383	else
2384	q_rico(k)=0.0
2385	endif
2386
2387	!T
2388	if(0 < zlay(k) .and. zlay(k) < 740) then
2389	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2390	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2391	t_rico(k)=292.0 + (278.0 - 292.0) / &
2392	& (4000 - 740) * (zlay(k) - 740)
2393	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2394	t_rico(k)=278.0 + (203.0 - 278.0) / &
2395	& (15000 - 4000) * (zlay(k) - 4000)
2396	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2397	t_rico(k)=203.0 + (194.0 - 203.0) / &
2398	& (17500 - 15000)* (zlay(k) - 15000)
2399	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2400	t_rico(k)=194.0 + (206.0 - 194.0) / &
2401	& (20000 - 17500)* (zlay(k) - 17500)
2402	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2403	t_rico(k)=206.0 + (270.0 - 206.0) / &
2404	& (60000 - 20000)* (zlay(k) - 20000)
2405	endif
2406
2407	! W
2408	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2409	w_rico(k)=- (0.005/2260) * zlay(k)
2410	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2411	w_rico(k)=- 0.005
2412	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2413	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2414	else
2415	w_rico(k)=0.0
2416	endif
2417
2418	! dThrz+dTsw0+dTlw0
2419	if(0 < zlay(k) .and. zlay(k) < 4000) then
2420	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ &
2421	& (864004000) zlay(k)
2422	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2423	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) / &
2424	& (86400(5000 - 4000)) (zlay(k) - 4000)
2425	else
2426	dth_dyn(k)=0.0
2427	endif
2428	! dQhrz
2429	if(0 < zlay(k) .and. zlay(k) < 3000) then
2430	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/ &
2431	& (864003000) (zlay(k))
2432	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2433	dqh_dyn(k)=0.345 / 86400
2434	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2435	dqh_dyn(k)=0.345 / 86400 + &
2436	& (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2437	else
2438	dqh_dyn(k)=0.0
2439	endif
2440
2441	!? if(play(k)>6e4) then
2442	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2443	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
2444	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2445	!? *(6e4-play(k))/(6e4-3e4)
2446	!? else
2447	!? ratqs0(1,k)=ratqshaut
2448	!? endif
2449
2450	enddo
2451
2452	do k=1,llm
2453	q_rico(k)=q_rico(k)/1e3
2454	dqh_dyn(k)=dqh_dyn(k)/1e3
2455	v_rico(k)=-3.8
2456	enddo
2457
2458	return
2459	end
2460
2461	!======================================================================
2462	SUBROUTINE interp_sandu_time(day,day1,annee_ref &
2463	& ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu &
2464	& ,nlev_sandu,ts_sandu,ts_prof)
2465	implicit none
2466
2467	!---------------------------------------------------------------------------------------
2468	! Time interpolation of a 2D field to the timestep corresponding to day
2469	!
2470	! day: current julian day (e.g. 717538.2)
2471	! day1: first day of the simulation
2472	! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
2473	! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
2474	!---------------------------------------------------------------------------------------
2475	! inputs:
2476	integer annee_ref
2477	integer nt_sandu,nlev_sandu
2478	integer year_ini_sandu
2479	real day, day1,day_ini_sandu,dt_sandu
2480	real ts_sandu(nt_sandu)
2481	! outputs:
2482	real ts_prof
2483	! local:
2484	integer it_sandu1, it_sandu2
2485	real timeit,time_sandu1,time_sandu2,frac
2486	! Check that initial day of the simulation consistent with SANDU period:
2487	if (annee_ref.ne.2006 ) then
2488	print*,'Pour SANDUREF, annee_ref doit etre 2006 '
2489	print*,'Changer annee_ref dans run.def'
2490	stop
2491	endif
2492	! if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
2493	! print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
2494	! print*,'Changer dayref dans run.def'
2495	! stop
2496	! endif
2497
2498	! Determine timestep relative to the 1st day of TOGA-COARE:
2499	! timeit=(day-day1)*86400.
2500	! if (annee_ref.eq.1992) then
2501	! timeit=(day-day_ini_sandu)*86400.
2502	! else
2503	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2504	! endif
2505	timeit=(day-day_ini_sandu)*86400
2506
2507	! Determine the closest observation times:
2508	it_sandu1=INT(timeit/dt_sandu)+1
2509	it_sandu2=it_sandu1 + 1
2510	time_sandu1=(it_sandu1-1)*dt_sandu
2511	time_sandu2=(it_sandu2-1)*dt_sandu
2512	print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
2513	print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2', &
2514	& it_sandu1,it_sandu2,time_sandu1,time_sandu2
2515
2516	if (it_sandu1 .ge. nt_sandu) then
2517	write(,) 'PB-stop: day, it_sandu1, it_sandu2, timeit: ' &
2518	& ,day,it_sandu1,it_sandu2,timeit/86400.
2519	stop
2520	endif
2521
2522	! time interpolation:
2523	frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
2524	frac=max(frac,0.0)
2525
2526	ts_prof = ts_sandu(it_sandu2) &
2527	& -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
2528
2529	print*, &
2530	&'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:', &
2531	&day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1, &
2532	&it_sandu2,ts_prof
2533
2534	return
2535	END
2536	!=====================================================================
2537	!-------------------------------------------------------------------------
2538	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu, &
2539	& sens,flat,adv_theta,rad_theta,adv_qt)
2540	implicit none
2541
2542	!-------------------------------------------------------------------------
2543	! Read ARM_CU case forcing data
2544	!-------------------------------------------------------------------------
2545
2546	integer nlev_armcu,nt_armcu
2547	real sens(nt_armcu),flat(nt_armcu)
2548	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2549	character*80 fich_armcu
2550
2551	integer ip
2552
2553	integer iy,im,id,ih,in
2554
2555	print*,'nlev_armcu',nlev_armcu
2556
2557	open(21,file=trim(fich_armcu),form='formatted')
2558	read(21,'(a)')
2559	do ip = 1, nt_armcu
2560	read(21,'(a)')
2561	read(21,'(a)')
2562	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip), &
2563	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
2564	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip), &
2565	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
2566	enddo
2567	close(21)
2568
2569	223 format(5i3,5f8.3)
2570
2571	return
2572	end
2573
2574	!=====================================================================
2575	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof &
2576	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
2577	& ,ht_prof,vt_prof,hq_prof,vq_prof &
2578	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
2579	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2580
2581	implicit none
2582
2583	#include "dimensions.h"
2584
2585	!-------------------------------------------------------------------------
2586	! Vertical interpolation of TOGA-COARE forcing data onto model levels
2587	!-------------------------------------------------------------------------
2588
2589	integer nlevmax
2590	parameter (nlevmax=41)
2591	integer nlev_toga,mxcalc
2592	! real play(llm), plev_prof(nlevmax)
2593	! real t_prof(nlevmax),q_prof(nlevmax)
2594	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2595	! real ht_prof(nlevmax),vt_prof(nlevmax)
2596	! real hq_prof(nlevmax),vq_prof(nlevmax)
2597
2598	real play(llm), plev_prof(nlev_toga)
2599	real t_prof(nlev_toga),q_prof(nlev_toga)
2600	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2601	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2602	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2603
2604	real t_mod(llm),q_mod(llm)
2605	real u_mod(llm),v_mod(llm), w_mod(llm)
2606	real ht_mod(llm),vt_mod(llm)
2607	real hq_mod(llm),vq_mod(llm)
2608
2609	integer l,k,k1,k2
2610	real frac,frac1,frac2,fact
2611
2612	do l = 1, llm
2613
2614	if (play(l).ge.plev_prof(nlev_toga)) then
2615
2616	mxcalc=l
2617	k1=0
2618	k2=0
2619
2620	if (play(l).le.plev_prof(1)) then
2621
2622	do k = 1, nlev_toga-1
2623	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2624	k1=k
2625	k2=k+1
2626	endif
2627	enddo
2628
2629	if (k1.eq.0 .or. k2.eq.0) then
2630	write(,) 'PB! k1, k2 = ',k1,k2
2631	write(,) 'l,play(l) = ',l,play(l)/100
2632	do k = 1, nlev_toga-1
2633	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2634	enddo
2635	endif
2636
2637	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2638	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2639	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2640	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2641	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2642	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2643	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2644	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2645	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2646	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2647
2648	else !play>plev_prof(1)
2649
2650	k1=1
2651	k2=2
2652	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2653	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2654	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2655	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
2656	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2657	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2658	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2659	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
2660	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
2661	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
2662	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
2663
2664	endif ! play.le.plev_prof(1)
2665
2666	else ! above max altitude of forcing file
2667
2668	!jyg
2669	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2670	fact = max(fact,0.) !jyg
2671	fact = exp(-fact) !jyg
2672	t_mod(l)= t_prof(nlev_toga) !jyg
2673	q_mod(l)= q_prof(nlev_toga)*fact !jyg
2674	u_mod(l)= u_prof(nlev_toga)*fact !jyg
2675	v_mod(l)= v_prof(nlev_toga)*fact !jyg
2676	w_mod(l)= 0.0 !jyg
2677	ht_mod(l)= ht_prof(nlev_toga) !jyg
2678	vt_mod(l)= vt_prof(nlev_toga) !jyg
2679	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
2680	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
2681
2682	endif ! play
2683
2684	enddo ! l
2685
2686	! do l = 1,llm
2687	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2688	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2689	! enddo
2690
2691	return
2692	end
2693
2694	!=====================================================================
2695	SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas &
2696	& ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &
2697	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
2698	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
2699	& ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &
2700	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
2701	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
2702
2703	implicit none
2704
2705	#include "dimensions.h"
2706
2707	!-------------------------------------------------------------------------
2708	! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
2709	!-------------------------------------------------------------------------
2710
2711	integer nlevmax
2712	parameter (nlevmax=41)
2713	integer nlev_cas,mxcalc
2714	! real play(llm), plev_prof(nlevmax)
2715	! real t_prof(nlevmax),q_prof(nlevmax)
2716	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2717	! real ht_prof(nlevmax),vt_prof(nlevmax)
2718	! real hq_prof(nlevmax),vq_prof(nlevmax)
2719
2720	real play(llm), plev_prof_cas(nlev_cas)
2721	real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
2722	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
2723	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
2724	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
2725	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
2726	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
2727	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
2728
2729	real t_mod_cas(llm),q_mod_cas(llm)
2730	real u_mod_cas(llm),v_mod_cas(llm)
2731	real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
2732	real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
2733	real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
2734	real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
2735	real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
2736
2737	integer l,k,k1,k2
2738	real frac,frac1,frac2,fact
2739
2740	do l = 1, llm
2741
2742	if (play(l).ge.plev_prof_cas(nlev_cas)) then
2743
2744	mxcalc=l
2745	k1=0
2746	k2=0
2747
2748	if (play(l).le.plev_prof_cas(1)) then
2749
2750	do k = 1, nlev_cas-1
2751	if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
2752	k1=k
2753	k2=k+1
2754	endif
2755	enddo
2756
2757	if (k1.eq.0 .or. k2.eq.0) then
2758	write(,) 'PB! k1, k2 = ',k1,k2
2759	write(,) 'l,play(l) = ',l,play(l)/100
2760	do k = 1, nlev_cas-1
2761	write(,) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
2762	enddo
2763	endif
2764
2765	frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
2766	t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
2767	q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
2768	u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
2769	v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
2770	ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
2771	vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
2772	w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
2773	du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
2774	hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
2775	vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
2776	dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
2777	hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
2778	vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
2779	dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
2780	ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
2781	vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
2782	dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
2783	hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
2784	vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
2785	dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
2786
2787	else !play>plev_prof_cas(1)
2788
2789	k1=1
2790	k2=2
2791	frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
2792	frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
2793	t_mod_cas(l)= frac1t_prof_cas(k1) - frac2t_prof_cas(k2)
2794	q_mod_cas(l)= frac1q_prof_cas(k1) - frac2q_prof_cas(k2)
2795	u_mod_cas(l)= frac1u_prof_cas(k1) - frac2u_prof_cas(k2)
2796	v_mod_cas(l)= frac1v_prof_cas(k1) - frac2v_prof_cas(k2)
2797	ug_mod_cas(l)= frac1ug_prof_cas(k1) - frac2ug_prof_cas(k2)
2798	vg_mod_cas(l)= frac1vg_prof_cas(k1) - frac2vg_prof_cas(k2)
2799	w_mod_cas(l)= frac1vitw_prof_cas(k1) - frac2vitw_prof_cas(k2)
2800	du_mod_cas(l)= frac1du_prof_cas(k1) - frac2du_prof_cas(k2)
2801	hu_mod_cas(l)= frac1hu_prof_cas(k1) - frac2hu_prof_cas(k2)
2802	vu_mod_cas(l)= frac1vu_prof_cas(k1) - frac2vu_prof_cas(k2)
2803	dv_mod_cas(l)= frac1dv_prof_cas(k1) - frac2dv_prof_cas(k2)
2804	hv_mod_cas(l)= frac1hv_prof_cas(k1) - frac2hv_prof_cas(k2)
2805	vv_mod_cas(l)= frac1vv_prof_cas(k1) - frac2vv_prof_cas(k2)
2806	dt_mod_cas(l)= frac1dt_prof_cas(k1) - frac2dt_prof_cas(k2)
2807	ht_mod_cas(l)= frac1ht_prof_cas(k1) - frac2ht_prof_cas(k2)
2808	vt_mod_cas(l)= frac1vt_prof_cas(k1) - frac2vt_prof_cas(k2)
2809	dq_mod_cas(l)= frac1dq_prof_cas(k1) - frac2dq_prof_cas(k2)
2810	hq_mod_cas(l)= frac1hq_prof_cas(k1) - frac2hq_prof_cas(k2)
2811	vq_mod_cas(l)= frac1vq_prof_cas(k1) - frac2vq_prof_cas(k2)
2812	dtrad_mod_cas(l)= frac1dtrad_prof_cas(k1) - frac2dtrad_prof_cas(k2)
2813
2814	endif ! play.le.plev_prof_cas(1)
2815
2816	else ! above max altitude of forcing file
2817
2818	!jyg
2819	fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
2820	fact = max(fact,0.) !jyg
2821	fact = exp(-fact) !jyg
2822	t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg
2823	q_mod_cas(l)= q_prof_cas(nlev_cas)*fact !jyg
2824	u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg
2825	v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg
2826	ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg
2827	vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg
2828	w_mod_cas(l)= 0.0 !jyg
2829	du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
2830	hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg
2831	vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg
2832	dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
2833	hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg
2834	vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg
2835	dt_mod_cas(l)= dt_prof_cas(nlev_cas)
2836	ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg
2837	vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg
2838	dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
2839	hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg
2840	vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg
2841	dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg
2842
2843	endif ! play
2844
2845	enddo ! l
2846
2847	! do l = 1,llm
2848	! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
2849	! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
2850	! enddo
2851
2852	return
2853	end
2854	!*****************************************************************************
2855	!=====================================================================
2856	SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof &
2857	& ,th_prof,qv_prof,u_prof,v_prof,o3_prof &
2858	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof &
2859	& ,th_mod,qv_mod,u_mod,v_mod,o3_mod &
2860	& ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
2861
2862	implicit none
2863
2864	#include "dimensions.h"
2865
2866	!-------------------------------------------------------------------------
2867	! Vertical interpolation of Dice forcing data onto model levels
2868	!-------------------------------------------------------------------------
2869
2870	integer nlevmax
2871	parameter (nlevmax=41)
2872	integer nlev_dice,mxcalc,nt_dice
2873
2874	real play(llm), plev_prof(nlev_dice)
2875	real th_prof(nlev_dice),qv_prof(nlev_dice)
2876	real u_prof(nlev_dice),v_prof(nlev_dice)
2877	real o3_prof(nlev_dice)
2878	real ht_prof(nlev_dice),hq_prof(nlev_dice)
2879	real hu_prof(nlev_dice),hv_prof(nlev_dice)
2880	real w_prof(nlev_dice),omega_prof(nlev_dice)
2881
2882	real th_mod(llm),qv_mod(llm)
2883	real u_mod(llm),v_mod(llm), o3_mod(llm)
2884	real ht_mod(llm),hq_mod(llm)
2885	real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
2886
2887	integer l,k,k1,k2,kp
2888	real aa,frac,frac1,frac2,fact
2889
2890	do l = 1, llm
2891
2892	if (play(l).ge.plev_prof(nlev_dice)) then
2893
2894	mxcalc=l
2895	k1=0
2896	k2=0
2897
2898	if (play(l).le.plev_prof(1)) then
2899
2900	do k = 1, nlev_dice-1
2901	if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
2902	k1=k
2903	k2=k+1
2904	endif
2905	enddo
2906
2907	if (k1.eq.0 .or. k2.eq.0) then
2908	write(,) 'PB! k1, k2 = ',k1,k2
2909	write(,) 'l,play(l) = ',l,play(l)/100
2910	do k = 1, nlev_dice-1
2911	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2912	enddo
2913	endif
2914
2915	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2916	th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
2917	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
2918	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2919	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2920	o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
2921	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2922	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2923	hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
2924	hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
2925	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2926	omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
2927
2928	else !play>plev_prof(1)
2929
2930	k1=1
2931	k2=2
2932	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2933	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2934	th_mod(l)= frac1th_prof(k1) - frac2th_prof(k2)
2935	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
2936	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2937	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2938	o3_mod(l)= frac1o3_prof(k1) - frac2o3_prof(k2)
2939	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
2940	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
2941	hu_mod(l)= frac1hu_prof(k1) - frac2hu_prof(k2)
2942	hv_mod(l)= frac1hv_prof(k1) - frac2hv_prof(k2)
2943	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2944	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
2945
2946	endif ! play.le.plev_prof(1)
2947
2948	else ! above max altitude of forcing file
2949
2950	!jyg
2951	fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
2952	fact = max(fact,0.) !jyg
2953	fact = exp(-fact) !jyg
2954	th_mod(l)= th_prof(nlev_dice) !jyg
2955	qv_mod(l)= qv_prof(nlev_dice)*fact !jyg
2956	u_mod(l)= u_prof(nlev_dice)*fact !jyg
2957	v_mod(l)= v_prof(nlev_dice)*fact !jyg
2958	o3_mod(l)= o3_prof(nlev_dice)*fact !jyg
2959	ht_mod(l)= ht_prof(nlev_dice) !jyg
2960	hq_mod(l)= hq_prof(nlev_dice)*fact !jyg
2961	hu_mod(l)= hu_prof(nlev_dice) !jyg
2962	hv_mod(l)= hv_prof(nlev_dice) !jyg
2963	w_mod(l)= 0. !jyg
2964	omega_mod(l)= 0. !jyg
2965
2966	endif ! play
2967
2968	enddo ! l
2969
2970	! do l = 1,llm
2971	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2972	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2973	! enddo
2974
2975	return
2976	end
2977
2978	!======================================================================
2979	SUBROUTINE interp_astex_time(day,day1,annee_ref &
2980	& ,year_ini_astex,day_ini_astex,nt_astex,dt_astex &
2981	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
2982	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
2983	& ,ufa_prof,vfa_prof)
2984	implicit none
2985
2986	!---------------------------------------------------------------------------------------
2987	! Time interpolation of a 2D field to the timestep corresponding to day
2988	!
2989	! day: current julian day (e.g. 717538.2)
2990	! day1: first day of the simulation
2991	! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
2992	! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
2993	!---------------------------------------------------------------------------------------
2994
2995	! inputs:
2996	integer annee_ref
2997	integer nt_astex,nlev_astex
2998	integer year_ini_astex
2999	real day, day1,day_ini_astex,dt_astex
3000	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
3001	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
3002	! outputs:
3003	real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
3004	! local:
3005	integer it_astex1, it_astex2
3006	real timeit,time_astex1,time_astex2,frac
3007
3008	! Check that initial day of the simulation consistent with ASTEX period:
3009	if (annee_ref.ne.1992 ) then
3010	print*,'Pour Astex, annee_ref doit etre 1992 '
3011	print*,'Changer annee_ref dans run.def'
3012	stop
3013	endif
3014	if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
3015	print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
3016	print*,'Changer dayref dans run.def'
3017	stop
3018	endif
3019
3020	! Determine timestep relative to the 1st day of TOGA-COARE:
3021	! timeit=(day-day1)*86400.
3022	! if (annee_ref.eq.1992) then
3023	! timeit=(day-day_ini_astex)*86400.
3024	! else
3025	! timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
3026	! endif
3027	timeit=(day-day_ini_astex)*86400
3028
3029	! Determine the closest observation times:
3030	it_astex1=INT(timeit/dt_astex)+1
3031	it_astex2=it_astex1 + 1
3032	time_astex1=(it_astex1-1)*dt_astex
3033	time_astex2=(it_astex2-1)*dt_astex
3034	print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
3035	print *,'it_astex1,it_astex2,time_astex1,time_astex2', &
3036	& it_astex1,it_astex2,time_astex1,time_astex2
3037
3038	if (it_astex1 .ge. nt_astex) then
3039	write(,) 'PB-stop: day, it_astex1, it_astex2, timeit: ' &
3040	& ,day,it_astex1,it_astex2,timeit/86400.
3041	stop
3042	endif
3043
3044	! time interpolation:
3045	frac=(time_astex2-timeit)/(time_astex2-time_astex1)
3046	frac=max(frac,0.0)
3047
3048	div_prof = div_astex(it_astex2) &
3049	& -frac*(div_astex(it_astex2)-div_astex(it_astex1))
3050	ts_prof = ts_astex(it_astex2) &
3051	& -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
3052	ug_prof = ug_astex(it_astex2) &
3053	& -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
3054	vg_prof = vg_astex(it_astex2) &
3055	& -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
3056	ufa_prof = ufa_astex(it_astex2) &
3057	& -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
3058	vfa_prof = vfa_astex(it_astex2) &
3059	& -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
3060
3061	print*, &
3062	&'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:', &
3063	&day,annee_ref,day_ini_astex,timeit/86400.,it_astex1, &
3064	&it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
3065
3066	return
3067	END
3068
3069	!======================================================================
3070	SUBROUTINE interp_toga_time(day,day1,annee_ref &
3071	& ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga &
3072	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
3073	& ,ht_toga,vt_toga,hq_toga,vq_toga &
3074	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
3075	& ,ht_prof,vt_prof,hq_prof,vq_prof)
3076	implicit none
3077
3078	!---------------------------------------------------------------------------------------
3079	! Time interpolation of a 2D field to the timestep corresponding to day
3080	!
3081	! day: current julian day (e.g. 717538.2)
3082	! day1: first day of the simulation
3083	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
3084	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
3085	!---------------------------------------------------------------------------------------
3086
3087	#include "compar1d.h"
3088
3089	! inputs:
3090	integer annee_ref
3091	integer nt_toga,nlev_toga
3092	integer year_ini_toga
3093	real day, day1,day_ini_toga,dt_toga
3094	real ts_toga(nt_toga)
3095	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
3096	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
3097	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
3098	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
3099	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
3100	! outputs:
3101	real ts_prof
3102	real plev_prof(nlev_toga),t_prof(nlev_toga)
3103	real q_prof(nlev_toga),u_prof(nlev_toga)
3104	real v_prof(nlev_toga),w_prof(nlev_toga)
3105	real ht_prof(nlev_toga),vt_prof(nlev_toga)
3106	real hq_prof(nlev_toga),vq_prof(nlev_toga)
3107	! local:
3108	integer it_toga1, it_toga2,k
3109	real timeit,time_toga1,time_toga2,frac
3110
3111
3112	if (forcing_type.eq.2) then
3113	! Check that initial day of the simulation consistent with TOGA-COARE period:
3114	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
3115	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
3116	print*,'Changer annee_ref dans run.def'
3117	stop
3118	endif
3119	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
3120	print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'
3121	print*,'Changer dayref dans run.def'
3122	stop
3123	endif
3124	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
3125	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
3126	print*,'Changer dayref ou nday dans run.def'
3127	stop
3128	endif
3129
3130	else if (forcing_type.eq.4) then
3131
3132	! Check that initial day of the simulation consistent with TWP-ICE period:
3133	if (annee_ref.ne.2006) then
3134	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
3135	print*,'Changer annee_ref dans run.def'
3136	stop
3137	endif
3138	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
3139	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
3140	print*,'Changer dayref dans run.def'
3141	stop
3142	endif
3143	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
3144	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
3145	print*,'Changer dayref ou nday dans run.def'
3146	stop
3147	endif
3148
3149	endif
3150
3151	! Determine timestep relative to the 1st day of TOGA-COARE:
3152	! timeit=(day-day1)*86400.
3153	! if (annee_ref.eq.1992) then
3154	! timeit=(day-day_ini_toga)*86400.
3155	! else
3156	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3157	! endif
3158	timeit=(day-day_ini_toga)*86400
3159
3160	! Determine the closest observation times:
3161	it_toga1=INT(timeit/dt_toga)+1
3162	it_toga2=it_toga1 + 1
3163	time_toga1=(it_toga1-1)*dt_toga
3164	time_toga2=(it_toga2-1)*dt_toga
3165
3166	if (it_toga1 .ge. nt_toga) then
3167	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: ' &
3168	& ,day,it_toga1,it_toga2,timeit/86400.
3169	stop
3170	endif
3171
3172	! time interpolation:
3173	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
3174	frac=max(frac,0.0)
3175
3176	ts_prof = ts_toga(it_toga2) &
3177	& -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
3178
3179	! print*,
3180	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
3181	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
3182
3183	do k=1,nlev_toga
3184	plev_prof(k) = 100.*(plev_toga(k,it_toga2) &
3185	& -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
3186	t_prof(k) = t_toga(k,it_toga2) &
3187	& -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
3188	q_prof(k) = q_toga(k,it_toga2) &
3189	& -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
3190	u_prof(k) = u_toga(k,it_toga2) &
3191	& -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
3192	v_prof(k) = v_toga(k,it_toga2) &
3193	& -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
3194	w_prof(k) = w_toga(k,it_toga2) &
3195	& -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
3196	ht_prof(k) = ht_toga(k,it_toga2) &
3197	& -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
3198	vt_prof(k) = vt_toga(k,it_toga2) &
3199	& -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
3200	hq_prof(k) = hq_toga(k,it_toga2) &
3201	& -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
3202	vq_prof(k) = vq_toga(k,it_toga2) &
3203	& -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
3204	enddo
3205
3206	return
3207	END
3208
3209	!======================================================================
3210	SUBROUTINE interp_dice_time(day,day1,annee_ref &
3211	& ,year_ini_dice,day_ini_dice,nt_dice,dt_dice &
3212	& ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice &
3213	& ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice &
3214	& ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
3215	& ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof &
3216	& ,ustar_prof,psurf_prof,ug_prof,vg_prof &
3217	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
3218	implicit none
3219
3220	!---------------------------------------------------------------------------------------
3221	! Time interpolation of a 2D field to the timestep corresponding to day
3222	!
3223	! day: current julian day (e.g. 717538.2)
3224	! day1: first day of the simulation
3225	! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
3226	! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
3227	!---------------------------------------------------------------------------------------
3228
3229	#include "compar1d.h"
3230
3231	! inputs:
3232	integer annee_ref
3233	integer nt_dice,nlev_dice
3234	integer year_ini_dice
3235	real day, day1,day_ini_dice,dt_dice
3236	real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
3237	real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
3238	real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
3239	real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
3240	real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
3241	real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
3242	! outputs:
3243	real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
3244	real ustar_prof,psurf_prof,ug_prof,vg_prof
3245	real ht_prof(nlev_dice),hq_prof(nlev_dice)
3246	real hu_prof(nlev_dice),hv_prof(nlev_dice)
3247	real w_prof(nlev_dice),omega_prof(nlev_dice)
3248	! local:
3249	integer it_dice1, it_dice2,k
3250	real timeit,time_dice1,time_dice2,frac
3251
3252
3253	if (forcing_type.eq.7) then
3254	! Check that initial day of the simulation consistent with Dice period:
3255	print *,'annee_ref=',annee_ref
3256	print *,'day1=',day1
3257	print *,'day_ini_dice=',day_ini_dice
3258	if (annee_ref.ne.1999) then
3259	print*,'Pour Dice, annee_ref doit etre 1999'
3260	print*,'Changer annee_ref dans run.def'
3261	stop
3262	endif
3263	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
3264	print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
3265	print*,'Changer dayref dans run.def',day1,day_ini_dice
3266	stop
3267	endif
3268	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
3269	print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
3270	print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
3271	stop
3272	endif
3273
3274	endif
3275
3276	! Determine timestep relative to the 1st day of TOGA-COARE:
3277	! timeit=(day-day1)*86400.
3278	! if (annee_ref.eq.1992) then
3279	! timeit=(day-day_ini_dice)*86400.
3280	! else
3281	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3282	! endif
3283	timeit=(day-day_ini_dice)*86400
3284
3285	! Determine the closest observation times:
3286	it_dice1=INT(timeit/dt_dice)+1
3287	it_dice2=it_dice1 + 1
3288	time_dice1=(it_dice1-1)*dt_dice
3289	time_dice2=(it_dice2-1)*dt_dice
3290
3291	if (it_dice1 .ge. nt_dice) then
3292	write(,) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
3293	stop
3294	endif
3295
3296	! time interpolation:
3297	frac=(time_dice2-timeit)/(time_dice2-time_dice1)
3298	frac=max(frac,0.0)
3299
3300	shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
3301	lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
3302	lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
3303	swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
3304	tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
3305	ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
3306	psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
3307	ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
3308	vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
3309
3310	! print*,
3311	! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
3312	! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
3313
3314	do k=1,nlev_dice
3315	ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
3316	hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
3317	hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
3318	hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
3319	w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
3320	omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
3321	enddo
3322
3323	return
3324	END
3325
3326	!======================================================================
3327	SUBROUTINE interp_gabls4_time(day,day1,annee_ref &
3328	& ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 &
3329	& ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 &
3330	& ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
3331	implicit none
3332
3333	!---------------------------------------------------------------------------------------
3334	! Time interpolation of a 2D field to the timestep corresponding to day
3335	!
3336	! day: current julian day
3337	! day1: first day of the simulation
3338	! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
3339	! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
3340	!---------------------------------------------------------------------------------------
3341
3342	#include "compar1d.h"
3343
3344	! inputs:
3345	integer annee_ref
3346	integer nt_gabls4,nlev_gabls4
3347	integer year_ini_gabls4
3348	real day, day1,day_ini_gabls4,dt_gabls4
3349	real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
3350	real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
3351	real tg_gabls4(nt_gabls4), tg_prof
3352	! outputs:
3353	real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
3354	real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
3355	! local:
3356	integer it_gabls41, it_gabls42,k
3357	real timeit,time_gabls41,time_gabls42,frac
3358
3359
3360
3361	! Check that initial day of the simulation consistent with gabls4 period:
3362	if (forcing_type.eq.8 ) then
3363	print *,'annee_ref=',annee_ref
3364	print *,'day1=',day1
3365	print *,'day_ini_gabls4=',day_ini_gabls4
3366	if (annee_ref.ne.2009) then
3367	print*,'Pour gabls4, annee_ref doit etre 2009'
3368	print*,'Changer annee_ref dans run.def'
3369	stop
3370	endif
3371	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
3372	print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
3373	print*,'Changer dayref dans run.def',day1,day_ini_gabls4
3374	stop
3375	endif
3376	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
3377	print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
3378	print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
3379	stop
3380	endif
3381	endif
3382
3383	timeit=(day-day_ini_gabls4)*86400
3384	print *,'day,day_ini_gabls4=',day,day_ini_gabls4
3385	print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
3386
3387	! Determine the closest observation times:
3388	it_gabls41=INT(timeit/dt_gabls4)+1
3389	it_gabls42=it_gabls41 + 1
3390	time_gabls41=(it_gabls41-1)*dt_gabls4
3391	time_gabls42=(it_gabls42-1)*dt_gabls4
3392
3393	if (it_gabls41 .ge. nt_gabls4) then
3394	write(,) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
3395	stop
3396	endif
3397
3398	! time interpolation:
3399	frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
3400	frac=max(frac,0.0)
3401
3402
3403	do k=1,nlev_gabls4
3404	ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
3405	vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
3406	ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
3407	hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
3408	enddo
3409	tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
3410	return
3411	END
3412
3413	!======================================================================
3414	SUBROUTINE interp_armcu_time(day,day1,annee_ref &
3415	& ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu &
3416	& ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu &
3417	& ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
3418	implicit none
3419
3420	!---------------------------------------------------------------------------------------
3421	! Time interpolation of a 2D field to the timestep corresponding to day
3422	!
3423	! day: current julian day (e.g. 717538.2)
3424	! day1: first day of the simulation
3425	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
3426	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
3427	! fs= sensible flux
3428	! fl= latent flux
3429	! at,rt,aqt= advective and radiative tendencies
3430	!---------------------------------------------------------------------------------------
3431
3432	! inputs:
3433	integer annee_ref
3434	integer nt_armcu,nlev_armcu
3435	integer year_ini_armcu
3436	real day, day1,day_ini_armcu,dt_armcu
3437	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
3438	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
3439	! outputs:
3440	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
3441	! local:
3442	integer it_armcu1, it_armcu2,k
3443	real timeit,time_armcu1,time_armcu2,frac
3444
3445	! Check that initial day of the simulation consistent with ARMCU period:
3446	if (annee_ref.ne.1997 ) then
3447	print*,'Pour ARMCU, annee_ref doit etre 1997 '
3448	print*,'Changer annee_ref dans run.def'
3449	stop
3450	endif
3451
3452	timeit=(day-day_ini_armcu)*86400
3453
3454	! Determine the closest observation times:
3455	it_armcu1=INT(timeit/dt_armcu)+1
3456	it_armcu2=it_armcu1 + 1
3457	time_armcu1=(it_armcu1-1)*dt_armcu
3458	time_armcu2=(it_armcu2-1)*dt_armcu
3459	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
3460	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', &
3461	& it_armcu1,it_armcu2,time_armcu1,time_armcu2
3462
3463	if (it_armcu1 .ge. nt_armcu) then
3464	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' &
3465	& ,day,it_armcu1,it_armcu2,timeit/86400.
3466	stop
3467	endif
3468
3469	! time interpolation:
3470	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
3471	frac=max(frac,0.0)
3472
3473	fs_prof = fs_armcu(it_armcu2) &
3474	& -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
3475	fl_prof = fl_armcu(it_armcu2) &
3476	& -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
3477	at_prof = at_armcu(it_armcu2) &
3478	& -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
3479	rt_prof = rt_armcu(it_armcu2) &
3480	& -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
3481	aqt_prof = aqt_armcu(it_armcu2) &
3482	& -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
3483
3484	print*, &
3485	&'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', &
3486	&day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, &
3487	&it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
3488
3489	return
3490	END
3491
3492	!=====================================================================
3493	subroutine readprofiles(nlev_max,kmax,ntrac,height, &
3494	& thlprof,qtprof,uprof, &
3495	& vprof,e12prof,ugprof,vgprof, &
3496	& wfls,dqtdxls,dqtdyls,dqtdtls, &
3497	& thlpcar,tracer,nt1,nt2)
3498	implicit none
3499
3500	integer nlev_max,kmax,kmax2,ntrac
3501	logical :: llesread = .true.
3502
3503	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), &
3504	& uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), &
3505	& ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), &
3506	& dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), &
3507	& thlpcar(nlev_max),tracer(nlev_max,ntrac)
3508
3509	real height1(nlev_max)
3510
3511	integer, parameter :: ilesfile=1
3512	integer :: ierr,k,itrac,nt1,nt2
3513
3514	if(.not.(llesread)) return
3515
3516	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3517	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3518	read (ilesfile,*) kmax
3519	do k=1,kmax
3520	read (ilesfile,*) height1(k),thlprof(k),qtprof (k), &
3521	& uprof (k),vprof (k),e12prof(k)
3522	enddo
3523	close(ilesfile)
3524
3525	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
3526	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
3527	read (ilesfile,*) kmax2
3528	if (kmax .ne. kmax2) then
3529	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3530	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3531	stop 'lecture profiles'
3532	endif
3533	do k=1,kmax
3534	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &
3535	& dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
3536	end do
3537	do k=1,kmax
3538	if (height(k) .ne. height1(k)) then
3539	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3540	print *, 'les niveaux different : ',k,height1(k), height(k)
3541	stop
3542	endif
3543	end do
3544	close(ilesfile)
3545
3546	open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
3547	if (ierr /= 0) then
3548	print*,'WARNING : trac.inp does not exist'
3549	else
3550	read (ilesfile,*) kmax2,nt1,nt2
3551	if (nt2>ntrac) then
3552	stop'Augmenter le nombre de traceurs dans traceur.def'
3553	endif
3554	if (kmax .ne. kmax2) then
3555	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3556	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3557	stop 'lecture profiles'
3558	endif
3559	do k=1,kmax
3560	read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
3561	end do
3562	close(ilesfile)
3563	endif
3564
3565	return
3566	end
3567	!======================================================================
3568	subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &
3569	& thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
3570	!======================================================================
3571	implicit none
3572
3573	integer nlev_max,kmax
3574	logical :: llesread = .true.
3575
3576	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3577	real thlprof(nlev_max)
3578	real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3579	real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
3580
3581	integer, parameter :: ilesfile=1
3582	integer :: k,ierr
3583
3584	if(.not.(llesread)) return
3585
3586	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3587	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3588	read (ilesfile,*) kmax
3589	do k=1,kmax
3590	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
3591	& qprof (k),uprof(k), vprof(k), wprof(k), &
3592	& omega (k),o3mmr(k)
3593	enddo
3594	close(ilesfile)
3595
3596	return
3597	end
3598
3599	!======================================================================
3600	subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &
3601	& thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
3602	!======================================================================
3603	implicit none
3604
3605	integer nlev_max,kmax
3606	logical :: llesread = .true.
3607
3608	real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &
3609	& thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &
3610	& qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &
3611	& wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
3612
3613	integer, parameter :: ilesfile=1
3614	integer :: ierr,k
3615
3616	if(.not.(llesread)) return
3617
3618	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3619	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3620	read (ilesfile,*) kmax
3621	do k=1,kmax
3622	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
3623	& qvprof (k),qlprof (k),qtprof (k), &
3624	& uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
3625	enddo
3626	close(ilesfile)
3627
3628	return
3629	end
3630
3631
3632
3633	!======================================================================
3634	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &
3635	& vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
3636	!======================================================================
3637	implicit none
3638
3639	integer nlev_max,kmax
3640	logical :: llesread = .true.
3641
3642	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3643	real thetaprof(nlev_max),rvprof(nlev_max)
3644	real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3645	real aprof(nlev_max+1),bprof(nlev_max+1)
3646
3647	integer, parameter :: ilesfile=1
3648	integer, parameter :: ifile=2
3649	integer :: ierr,jtot,k
3650
3651	if(.not.(llesread)) return
3652
3653	! Read profiles at full levels
3654	IF(nlev_max.EQ.19) THEN
3655	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
3656	print *,'On ouvre prof.inp.19'
3657	ELSE
3658	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
3659	print *,'On ouvre prof.inp.40'
3660	ENDIF
3661	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3662	read (ilesfile,*) kmax
3663	do k=1,kmax
3664	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &
3665	& thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
3666	enddo
3667	close(ilesfile)
3668
3669	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
3670	IF(nlev_max.EQ.19) THEN
3671	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
3672	print *,'On ouvre proh.inp.19'
3673	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
3674	ELSE
3675	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
3676	print *,'On ouvre proh.inp.40'
3677	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
3678	ENDIF
3679	read (ifile,*) kmax
3680	do k=1,kmax
3681	read (ifile,*) jtot,aprof(k),bprof(k)
3682	enddo
3683	close(ifile)
3684
3685	return
3686	end
3687
3688	!=====================================================================
3689	subroutine read_fire(fich_fire,nlevel,ntime &
3690	& ,zz,thl,qt,u,v,tke &
3691	& ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3692
3693	!program reading forcings of the FIRE case study
3694
3695
3696	implicit none
3697
3698	#include "netcdf.inc"
3699
3700	integer ntime,nlevel
3701	character*80 :: fich_fire
3702	real*8 zz(nlevel)
3703
3704	real*8 thl(nlevel)
3705	real*8 qt(nlevel),u(nlevel)
3706	real*8 v(nlevel),tke(nlevel)
3707	real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3708	real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3709	real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3710
3711	integer nid, ierr
3712	integer nbvar3d
3713	parameter(nbvar3d=30)
3714	integer var3didin(nbvar3d)
3715
3716	ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3717	if (ierr.NE.NF_NOERR) then
3718	write(,) 'ERROR: Pb opening forcings nc file '
3719	write(,) NF_STRERROR(ierr)
3720	stop ""
3721	endif
3722
3723
3724	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
3725	if(ierr/=NF_NOERR) then
3726	write(,) NF_STRERROR(ierr)
3727	stop 'lev'
3728	endif
3729
3730
3731	ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3732	if(ierr/=NF_NOERR) then
3733	write(,) NF_STRERROR(ierr)
3734	stop 'temp'
3735	endif
3736
3737	ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3738	if(ierr/=NF_NOERR) then
3739	write(,) NF_STRERROR(ierr)
3740	stop 'qv'
3741	endif
3742
3743	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3744	if(ierr/=NF_NOERR) then
3745	write(,) NF_STRERROR(ierr)
3746	stop 'u'
3747	endif
3748
3749	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3750	if(ierr/=NF_NOERR) then
3751	write(,) NF_STRERROR(ierr)
3752	stop 'v'
3753	endif
3754
3755	ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3756	if(ierr/=NF_NOERR) then
3757	write(,) NF_STRERROR(ierr)
3758	stop 'tke'
3759	endif
3760
3761	ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3762	if(ierr/=NF_NOERR) then
3763	write(,) NF_STRERROR(ierr)
3764	stop 'ug'
3765	endif
3766
3767	ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3768	if(ierr/=NF_NOERR) then
3769	write(,) NF_STRERROR(ierr)
3770	stop 'vg'
3771	endif
3772
3773	ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3774	if(ierr/=NF_NOERR) then
3775	write(,) NF_STRERROR(ierr)
3776	stop 'wls'
3777	endif
3778
3779	ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3780	if(ierr/=NF_NOERR) then
3781	write(,) NF_STRERROR(ierr)
3782	stop 'dqtdx'
3783	endif
3784
3785	ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3786	if(ierr/=NF_NOERR) then
3787	write(,) NF_STRERROR(ierr)
3788	stop 'dqtdy'
3789	endif
3790
3791	ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3792	if(ierr/=NF_NOERR) then
3793	write(,) NF_STRERROR(ierr)
3794	stop 'dqtdt'
3795	endif
3796
3797	ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3798	if(ierr/=NF_NOERR) then
3799	write(,) NF_STRERROR(ierr)
3800	stop 'thl_rad'
3801	endif
3802	!dimensions lecture
3803	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3804
3805	#ifdef NC_DOUBLE
3806	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3807	#else
3808	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3809	#endif
3810	if(ierr/=NF_NOERR) then
3811	write(,) NF_STRERROR(ierr)
3812	stop "getvarup"
3813	endif
3814	! write(,)'lecture z ok',zz
3815
3816	#ifdef NC_DOUBLE
3817	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3818	#else
3819	ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3820	#endif
3821	if(ierr/=NF_NOERR) then
3822	write(,) NF_STRERROR(ierr)
3823	stop "getvarup"
3824	endif
3825	! write(,)'lecture thl ok',thl
3826
3827	#ifdef NC_DOUBLE
3828	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3829	#else
3830	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3831	#endif
3832	if(ierr/=NF_NOERR) then
3833	write(,) NF_STRERROR(ierr)
3834	stop "getvarup"
3835	endif
3836	! write(,)'lecture qt ok',qt
3837
3838	#ifdef NC_DOUBLE
3839	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3840	#else
3841	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3842	#endif
3843	if(ierr/=NF_NOERR) then
3844	write(,) NF_STRERROR(ierr)
3845	stop "getvarup"
3846	endif
3847	! write(,)'lecture u ok',u
3848
3849	#ifdef NC_DOUBLE
3850	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3851	#else
3852	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3853	#endif
3854	if(ierr/=NF_NOERR) then
3855	write(,) NF_STRERROR(ierr)
3856	stop "getvarup"
3857	endif
3858	! write(,)'lecture v ok',v
3859
3860	#ifdef NC_DOUBLE
3861	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3862	#else
3863	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3864	#endif
3865	if(ierr/=NF_NOERR) then
3866	write(,) NF_STRERROR(ierr)
3867	stop "getvarup"
3868	endif
3869	! write(,)'lecture tke ok',tke
3870
3871	#ifdef NC_DOUBLE
3872	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3873	#else
3874	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3875	#endif
3876	if(ierr/=NF_NOERR) then
3877	write(,) NF_STRERROR(ierr)
3878	stop "getvarup"
3879	endif
3880	! write(,)'lecture ug ok',ug
3881
3882	#ifdef NC_DOUBLE
3883	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3884	#else
3885	ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3886	#endif
3887	if(ierr/=NF_NOERR) then
3888	write(,) NF_STRERROR(ierr)
3889	stop "getvarup"
3890	endif
3891	! write(,)'lecture vg ok',vg
3892
3893	#ifdef NC_DOUBLE
3894	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3895	#else
3896	ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3897	#endif
3898	if(ierr/=NF_NOERR) then
3899	write(,) NF_STRERROR(ierr)
3900	stop "getvarup"
3901	endif
3902	! write(,)'lecture wls ok',wls
3903
3904	#ifdef NC_DOUBLE
3905	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3906	#else
3907	ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3908	#endif
3909	if(ierr/=NF_NOERR) then
3910	write(,) NF_STRERROR(ierr)
3911	stop "getvarup"
3912	endif
3913	! write(,)'lecture dqtdx ok',dqtdx
3914
3915	#ifdef NC_DOUBLE
3916	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3917	#else
3918	ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3919	#endif
3920	if(ierr/=NF_NOERR) then
3921	write(,) NF_STRERROR(ierr)
3922	stop "getvarup"
3923	endif
3924	! write(,)'lecture dqtdy ok',dqtdy
3925
3926	#ifdef NC_DOUBLE
3927	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3928	#else
3929	ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3930	#endif
3931	if(ierr/=NF_NOERR) then
3932	write(,) NF_STRERROR(ierr)
3933	stop "getvarup"
3934	endif
3935	! write(,)'lecture dqtdt ok',dqtdt
3936
3937	#ifdef NC_DOUBLE
3938	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3939	#else
3940	ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3941	#endif
3942	if(ierr/=NF_NOERR) then
3943	write(,) NF_STRERROR(ierr)
3944	stop "getvarup"
3945	endif
3946	! write(,)'lecture thl_rad ok',thl_rad
3947
3948	return
3949	end subroutine read_fire
3950	!=====================================================================
3951	subroutine read_dice(fich_dice,nlevel,ntime &
3952	& ,zz,pres,t,qv,u,v,o3 &
3953	& ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg &
3954	& ,hadvt,hadvq,hadvu,hadvv,w,omega)
3955
3956	!program reading initial profils and forcings of the Dice case study
3957
3958
3959	implicit none
3960
3961	#include "netcdf.inc"
3962	#include "YOMCST.h"
3963
3964	integer ntime,nlevel
3965	integer l,k
3966	character*80 :: fich_dice
3967	real*8 time(ntime)
3968	real*8 zz(nlevel)
3969
3970	real*8 th(nlevel),pres(nlevel),t(nlevel)
3971	real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
3972	real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
3973	real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
3974	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
3975	real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
3976	real*8 pzero
3977
3978	integer nid, ierr
3979	integer nbvar3d
3980	parameter(nbvar3d=30)
3981	integer var3didin(nbvar3d)
3982
3983	pzero=100000.
3984	ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
3985	if (ierr.NE.NF_NOERR) then
3986	write(,) 'ERROR: Pb opening forcings nc file '
3987	write(,) NF_STRERROR(ierr)
3988	stop ""
3989	endif
3990
3991
3992	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
3993	if(ierr/=NF_NOERR) then
3994	write(,) NF_STRERROR(ierr)
3995	stop 'height'
3996	endif
3997
3998	ierr=NF_INQ_VARID(nid,"pf",var3didin(11))
3999	if(ierr/=NF_NOERR) then
4000	write(,) NF_STRERROR(ierr)
4001	stop 'pf'
4002	endif
4003
4004	ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
4005	if(ierr/=NF_NOERR) then
4006	write(,) NF_STRERROR(ierr)
4007	stop 'theta'
4008	endif
4009
4010	ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
4011	if(ierr/=NF_NOERR) then
4012	write(,) NF_STRERROR(ierr)
4013	stop 'qv'
4014	endif
4015
4016	ierr=NF_INQ_VARID(nid,"u",var3didin(14))
4017	if(ierr/=NF_NOERR) then
4018	write(,) NF_STRERROR(ierr)
4019	stop 'u'
4020	endif
4021
4022	ierr=NF_INQ_VARID(nid,"v",var3didin(15))
4023	if(ierr/=NF_NOERR) then
4024	write(,) NF_STRERROR(ierr)
4025	stop 'v'
4026	endif
4027
4028	ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
4029	if(ierr/=NF_NOERR) then
4030	write(,) NF_STRERROR(ierr)
4031	stop 'o3'
4032	endif
4033
4034	ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
4035	if(ierr/=NF_NOERR) then
4036	write(,) NF_STRERROR(ierr)
4037	stop 'shf'
4038	endif
4039
4040	ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
4041	if(ierr/=NF_NOERR) then
4042	write(,) NF_STRERROR(ierr)
4043	stop 'lhf'
4044	endif
4045
4046	ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
4047	if(ierr/=NF_NOERR) then
4048	write(,) NF_STRERROR(ierr)
4049	stop 'lwup'
4050	endif
4051
4052	ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
4053	if(ierr/=NF_NOERR) then
4054	write(,) NF_STRERROR(ierr)
4055	stop 'dqtdx'
4056	endif
4057
4058	ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
4059	if(ierr/=NF_NOERR) then
4060	write(,) NF_STRERROR(ierr)
4061	stop 'Tg'
4062	endif
4063
4064	ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
4065	if(ierr/=NF_NOERR) then
4066	write(,) NF_STRERROR(ierr)
4067	stop 'ustar'
4068	endif
4069
4070	ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
4071	if(ierr/=NF_NOERR) then
4072	write(,) NF_STRERROR(ierr)
4073	stop 'psurf'
4074	endif
4075
4076	ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
4077	if(ierr/=NF_NOERR) then
4078	write(,) NF_STRERROR(ierr)
4079	stop 'Ug'
4080	endif
4081
4082	ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
4083	if(ierr/=NF_NOERR) then
4084	write(,) NF_STRERROR(ierr)
4085	stop 'Vg'
4086	endif
4087
4088	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
4089	if(ierr/=NF_NOERR) then
4090	write(,) NF_STRERROR(ierr)
4091	stop 'hadvT'
4092	endif
4093
4094	ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
4095	if(ierr/=NF_NOERR) then
4096	write(,) NF_STRERROR(ierr)
4097	stop 'hadvq'
4098	endif
4099
4100	ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
4101	if(ierr/=NF_NOERR) then
4102	write(,) NF_STRERROR(ierr)
4103	stop 'hadvu'
4104	endif
4105
4106	ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
4107	if(ierr/=NF_NOERR) then
4108	write(,) NF_STRERROR(ierr)
4109	stop 'hadvv'
4110	endif
4111
4112	ierr=NF_INQ_VARID(nid,"w",var3didin(21))
4113	if(ierr/=NF_NOERR) then
4114	write(,) NF_STRERROR(ierr)
4115	stop 'w'
4116	endif
4117
4118	ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
4119	if(ierr/=NF_NOERR) then
4120	write(,) NF_STRERROR(ierr)
4121	stop 'omega'
4122	endif
4123	!dimensions lecture
4124	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
4125
4126	#ifdef NC_DOUBLE
4127	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
4128	#else
4129	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
4130	#endif
4131	if(ierr/=NF_NOERR) then
4132	write(,) NF_STRERROR(ierr)
4133	stop "getvarup"
4134	endif
4135	! write(,)'lecture zz ok',zz
4136
4137	#ifdef NC_DOUBLE
4138	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
4139	#else
4140	ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
4141	#endif
4142	if(ierr/=NF_NOERR) then
4143	write(,) NF_STRERROR(ierr)
4144	stop "getvarup"
4145	endif
4146	! write(,)'lecture pres ok',pres
4147
4148	#ifdef NC_DOUBLE
4149	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
4150	#else
4151	ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
4152	#endif
4153	if(ierr/=NF_NOERR) then
4154	write(,) NF_STRERROR(ierr)
4155	stop "getvarup"
4156	endif
4157	! write(,)'lecture th ok',th
4158	do k=1,nlevel
4159	t(k)=th(k)(pres(k)/pzero)*rkappa
4160	enddo
4161
4162	#ifdef NC_DOUBLE
4163	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
4164	#else
4165	ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
4166	#endif
4167	if(ierr/=NF_NOERR) then
4168	write(,) NF_STRERROR(ierr)
4169	stop "getvarup"
4170	endif
4171	! write(,)'lecture qv ok',qv
4172
4173	#ifdef NC_DOUBLE
4174	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
4175	#else
4176	ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
4177	#endif
4178	if(ierr/=NF_NOERR) then
4179	write(,) NF_STRERROR(ierr)
4180	stop "getvarup"
4181	endif
4182	! write(,)'lecture u ok',u
4183
4184	#ifdef NC_DOUBLE
4185	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
4186	#else
4187	ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
4188	#endif
4189	if(ierr/=NF_NOERR) then
4190	write(,) NF_STRERROR(ierr)
4191	stop "getvarup"
4192	endif
4193	! write(,)'lecture v ok',v
4194
4195	#ifdef NC_DOUBLE
4196	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
4197	#else
4198	ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
4199	#endif
4200	if(ierr/=NF_NOERR) then
4201	write(,) NF_STRERROR(ierr)
4202	stop "getvarup"
4203	endif
4204	! write(,)'lecture o3 ok',o3
4205
4206	#ifdef NC_DOUBLE
4207	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
4208	#else
4209	ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
4210	#endif
4211	if(ierr/=NF_NOERR) then
4212	write(,) NF_STRERROR(ierr)
4213	stop "getvarup"
4214	endif
4215	! write(,)'lecture shf ok',shf
4216
4217	#ifdef NC_DOUBLE
4218	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
4219	#else
4220	ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
4221	#endif
4222	if(ierr/=NF_NOERR) then
4223	write(,) NF_STRERROR(ierr)
4224	stop "getvarup"
4225	endif
4226	! write(,)'lecture lhf ok',lhf
4227
4228	#ifdef NC_DOUBLE
4229	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
4230	#else
4231	ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
4232	#endif
4233	if(ierr/=NF_NOERR) then
4234	write(,) NF_STRERROR(ierr)
4235	stop "getvarup"
4236	endif
4237	! write(,)'lecture lwup ok',lwup
4238
4239	#ifdef NC_DOUBLE
4240	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
4241	#else
4242	ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
4243	#endif
4244	if(ierr/=NF_NOERR) then
4245	write(,) NF_STRERROR(ierr)
4246	stop "getvarup"
4247	endif
4248	! write(,)'lecture swup ok',swup
4249
4250	#ifdef NC_DOUBLE
4251	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
4252	#else
4253	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
4254	#endif
4255	if(ierr/=NF_NOERR) then
4256	write(,) NF_STRERROR(ierr)
4257	stop "getvarup"
4258	endif
4259	! write(,)'lecture tg ok',tg
4260
4261	#ifdef NC_DOUBLE
4262	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
4263	#else
4264	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
4265	#endif
4266	if(ierr/=NF_NOERR) then
4267	write(,) NF_STRERROR(ierr)
4268	stop "getvarup"
4269	endif
4270	! write(,)'lecture ustar ok',ustar
4271
4272	#ifdef NC_DOUBLE
4273	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
4274	#else
4275	ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
4276	#endif
4277	if(ierr/=NF_NOERR) then
4278	write(,) NF_STRERROR(ierr)
4279	stop "getvarup"
4280	endif
4281	! write(,)'lecture psurf ok',psurf
4282
4283	#ifdef NC_DOUBLE
4284	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
4285	#else
4286	ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
4287	#endif
4288	if(ierr/=NF_NOERR) then
4289	write(,) NF_STRERROR(ierr)
4290	stop "getvarup"
4291	endif
4292	! write(,)'lecture ug ok',ug
4293
4294	#ifdef NC_DOUBLE
4295	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
4296	#else
4297	ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
4298	#endif
4299	if(ierr/=NF_NOERR) then
4300	write(,) NF_STRERROR(ierr)
4301	stop "getvarup"
4302	endif
4303	! write(,)'lecture vg ok',vg
4304
4305	#ifdef NC_DOUBLE
4306	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
4307	#else
4308	ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
4309	#endif
4310	if(ierr/=NF_NOERR) then
4311	write(,) NF_STRERROR(ierr)
4312	stop "getvarup"
4313	endif
4314	! write(,)'lecture hadvt ok',hadvt
4315
4316	#ifdef NC_DOUBLE
4317	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
4318	#else
4319	ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
4320	#endif
4321	if(ierr/=NF_NOERR) then
4322	write(,) NF_STRERROR(ierr)
4323	stop "getvarup"
4324	endif
4325	! write(,)'lecture hadvq ok',hadvq
4326
4327	#ifdef NC_DOUBLE
4328	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
4329	#else
4330	ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
4331	#endif
4332	if(ierr/=NF_NOERR) then
4333	write(,) NF_STRERROR(ierr)
4334	stop "getvarup"
4335	endif
4336	! write(,)'lecture hadvu ok',hadvu
4337
4338	#ifdef NC_DOUBLE
4339	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
4340	#else
4341	ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
4342	#endif
4343	if(ierr/=NF_NOERR) then
4344	write(,) NF_STRERROR(ierr)
4345	stop "getvarup"
4346	endif
4347	! write(,)'lecture hadvv ok',hadvv
4348
4349	#ifdef NC_DOUBLE
4350	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
4351	#else
4352	ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
4353	#endif
4354	if(ierr/=NF_NOERR) then
4355	write(,) NF_STRERROR(ierr)
4356	stop "getvarup"
4357	endif
4358	! write(,)'lecture w ok',w
4359
4360	#ifdef NC_DOUBLE
4361	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
4362	#else
4363	ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
4364	#endif
4365	if(ierr/=NF_NOERR) then
4366	write(,) NF_STRERROR(ierr)
4367	stop "getvarup"
4368	endif
4369	! write(,)'lecture omega ok',omega
4370
4371	return
4372	end subroutine read_dice
4373	!=====================================================================
4374	subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol &
4375	& ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
4376
4377	!program reading initial profils and forcings of the Gabls4 case study
4378
4379
4380	implicit none
4381
4382	#include "netcdf.inc"
4383
4384	integer ntime,nlevel,nsol
4385	integer l,k
4386	character*80 :: fich_gabls4
4387	real*8 time(ntime)
4388
4389	! ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees
4390	! dans un ordre inverse par rapport a la convention LMDZ
4391	! ==> il faut tout inverser (MPL 20141024)
4392	! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
4393	real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
4394	real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
4395	real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
4396
4397	real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
4398	real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
4399	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
4400
4401	real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
4402	real*8 tg(ntime)
4403	integer nid, ierr
4404	integer nbvar3d
4405	parameter(nbvar3d=30)
4406	integer var3didin(nbvar3d)
4407
4408	ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
4409	if (ierr.NE.NF_NOERR) then
4410	write(,) 'ERROR: Pb opening forcings nc file '
4411	write(,) NF_STRERROR(ierr)
4412	stop ""
4413	endif
4414
4415
4416	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
4417	if(ierr/=NF_NOERR) then
4418	write(,) NF_STRERROR(ierr)
4419	stop 'height'
4420	endif
4421
4422	ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
4423	if(ierr/=NF_NOERR) then
4424	write(,) NF_STRERROR(ierr)
4425	stop 'depth_sn'
4426	endif
4427
4428	ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
4429	if(ierr/=NF_NOERR) then
4430	write(,) NF_STRERROR(ierr)
4431	stop 'Ug'
4432	endif
4433
4434	ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
4435	if(ierr/=NF_NOERR) then
4436	write(,) NF_STRERROR(ierr)
4437	stop 'Vg'
4438	endif
4439	ierr=NF_INQ_VARID(nid,"pf",var3didin(5))
4440	if(ierr/=NF_NOERR) then
4441	write(,) NF_STRERROR(ierr)
4442	stop 'pf'
4443	endif
4444
4445	ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
4446	if(ierr/=NF_NOERR) then
4447	write(,) NF_STRERROR(ierr)
4448	stop 'theta'
4449	endif
4450
4451	ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
4452	if(ierr/=NF_NOERR) then
4453	write(,) NF_STRERROR(ierr)
4454	stop 'tempe'
4455	endif
4456
4457	ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
4458	if(ierr/=NF_NOERR) then
4459	write(,) NF_STRERROR(ierr)
4460	stop 'qv'
4461	endif
4462
4463	ierr=NF_INQ_VARID(nid,"u",var3didin(9))
4464	if(ierr/=NF_NOERR) then
4465	write(,) NF_STRERROR(ierr)
4466	stop 'u'
4467	endif
4468
4469	ierr=NF_INQ_VARID(nid,"v",var3didin(10))
4470	if(ierr/=NF_NOERR) then
4471	write(,) NF_STRERROR(ierr)
4472	stop 'v'
4473	endif
4474
4475	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
4476	if(ierr/=NF_NOERR) then
4477	write(,) NF_STRERROR(ierr)
4478	stop 'hadvt'
4479	endif
4480
4481	ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
4482	if(ierr/=NF_NOERR) then
4483	write(,) NF_STRERROR(ierr)
4484	stop 'hadvq'
4485	endif
4486
4487	ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
4488	if(ierr/=NF_NOERR) then
4489	write(,) NF_STRERROR(ierr)
4490	stop 'tsnow'
4491	endif
4492
4493	ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
4494	if(ierr/=NF_NOERR) then
4495	write(,) NF_STRERROR(ierr)
4496	stop 'snow_density'
4497	endif
4498
4499	ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
4500	if(ierr/=NF_NOERR) then
4501	write(,) NF_STRERROR(ierr)
4502	stop 'Tg'
4503	endif
4504
4505
4506	!dimensions lecture
4507	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
4508
4509	#ifdef NC_DOUBLE
4510	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)
4511	#else
4512	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)
4513	#endif
4514	if(ierr/=NF_NOERR) then
4515	write(,) NF_STRERROR(ierr)
4516	stop "getvarup"
4517	endif
4518
4519	#ifdef NC_DOUBLE
4520	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)
4521	#else
4522	ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)
4523	#endif
4524	if(ierr/=NF_NOERR) then
4525	write(,) NF_STRERROR(ierr)
4526	stop "getvarup"
4527	endif
4528
4529	#ifdef NC_DOUBLE
4530	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)
4531	#else
4532	ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)
4533	#endif
4534	if(ierr/=NF_NOERR) then
4535	write(,) NF_STRERROR(ierr)
4536	stop "getvarup"
4537	endif
4538
4539	#ifdef NC_DOUBLE
4540	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)
4541	#else
4542	ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)
4543	#endif
4544	if(ierr/=NF_NOERR) then
4545	write(,) NF_STRERROR(ierr)
4546	stop "getvarup"
4547	endif
4548
4549	#ifdef NC_DOUBLE
4550	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)
4551	#else
4552	ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)
4553	#endif
4554	if(ierr/=NF_NOERR) then
4555	write(,) NF_STRERROR(ierr)
4556	stop "getvarup"
4557	endif
4558
4559	#ifdef NC_DOUBLE
4560	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)
4561	#else
4562	ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)
4563	#endif
4564	if(ierr/=NF_NOERR) then
4565	write(,) NF_STRERROR(ierr)
4566	stop "getvarup"
4567	endif
4568
4569	#ifdef NC_DOUBLE
4570	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)
4571	#else
4572	ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)
4573	#endif
4574	if(ierr/=NF_NOERR) then
4575	write(,) NF_STRERROR(ierr)
4576	stop "getvarup"
4577	endif
4578
4579	#ifdef NC_DOUBLE
4580	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)
4581	#else
4582	ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)
4583	#endif
4584	if(ierr/=NF_NOERR) then
4585	write(,) NF_STRERROR(ierr)
4586	stop "getvarup"
4587	endif
4588
4589	#ifdef NC_DOUBLE
4590	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)
4591	#else
4592	ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)
4593	#endif
4594	if(ierr/=NF_NOERR) then
4595	write(,) NF_STRERROR(ierr)
4596	stop "getvarup"
4597	endif
4598
4599	#ifdef NC_DOUBLE
4600	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)
4601	#else
4602	ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)
4603	#endif
4604	if(ierr/=NF_NOERR) then
4605	write(,) NF_STRERROR(ierr)
4606	stop "getvarup"
4607	endif
4608
4609	#ifdef NC_DOUBLE
4610	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)
4611	#else
4612	ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)
4613	#endif
4614	if(ierr/=NF_NOERR) then
4615	write(,) NF_STRERROR(ierr)
4616	stop "getvarup"
4617	endif
4618
4619	#ifdef NC_DOUBLE
4620	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)
4621	#else
4622	ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)
4623	#endif
4624	if(ierr/=NF_NOERR) then
4625	write(,) NF_STRERROR(ierr)
4626	stop "getvarup"
4627	endif
4628
4629	#ifdef NC_DOUBLE
4630	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)
4631	#else
4632	ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)
4633	#endif
4634	if(ierr/=NF_NOERR) then
4635	write(,) NF_STRERROR(ierr)
4636	stop "getvarup"
4637	endif
4638
4639	#ifdef NC_DOUBLE
4640	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)
4641	#else
4642	ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)
4643	#endif
4644	if(ierr/=NF_NOERR) then
4645	write(,) NF_STRERROR(ierr)
4646	stop "getvarup"
4647	endif
4648
4649	#ifdef NC_DOUBLE
4650	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)
4651	#else
4652	ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)
4653	#endif
4654	if(ierr/=NF_NOERR) then
4655	write(,) NF_STRERROR(ierr)
4656	stop "getvarup"
4657	endif
4658
4659	! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
4660	do k=1,nlevel
4661	zz(k)=zz_i(nlevel+1-k)
4662	ug(k,:)=ug_i(nlevel+1-k,:)
4663	vg(k,:)=vg_i(nlevel+1-k,:)
4664	pf(k)=pf_i(nlevel+1-k)
4665	print *,'pf=',pf(k)
4666	th(k)=th_i(nlevel+1-k)
4667	t(k)=t_i(nlevel+1-k)
4668	qv(k)=qv_i(nlevel+1-k)
4669	u(k)=u_i(nlevel+1-k)
4670	v(k)=v_i(nlevel+1-k)
4671	hadvt(k,:)=hadvt_i(nlevel+1-k,:)
4672	hadvq(k,:)=hadvq_i(nlevel+1-k,:)
4673	enddo
4674	return
4675	end subroutine read_gabls4
4676	!=====================================================================
4677
4678	! Reads CIRC input files
4679
4680	SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
4681
4682	parameter (ncm_1=49180)
4683	#include "YOMCST.h"
4684
4685	real albsfc(ncm_1), albsfc_w(ncm_1)
4686	real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
4687	reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
4688	real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
4689	real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
4690	real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
4691	real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
4692	o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
4693	! za= zenital angle
4694	! sza= cosinus angle zenital
4695	real wavn(ncm_1), ssf(ncm_1),za,sza
4696	integer nlev
4697
4698
4699	! Open the files
4700
4701	open (11, file='Tsfc_sza_nlev_case.txt', status='old')
4702	open (12, file='level_input_case.txt', status='old')
4703	open (13, file='layer_input_case.txt', status='old')
4704	open (14, file='aerosol_input_case.txt', status='old')
4705	open (15, file='cloud_input_case.txt', status='old')
4706	open (16, file='sfcalbedo_input_case.txt', status='old')
4707
4708	! Read scalar information
4709	do iskip=1,5
4710	read (11, *)
4711	enddo
4712	read (11, '(i8)') nlev
4713	read (11, '(f10.2)') tsfc
4714	read (11, '(f10.2)') za
4715	read (11, '(f10.4)') sw_dn_toa
4716	sza=cos(za/180.*RPI)
4717	print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
4718	close(11)
4719
4720	! Read level information
4721	read (12, *)
4722	do il=1,nlev
4723	read (12, 302) ilev, z(il), p(il), t(il)
4724	z(il)=z(il)*1000. ! z donne en km
4725	p(il)=p(il)*100. ! p donne en mb
4726	enddo
4727	302 format (i8, f8.3, 2f9.2)
4728	close(12)
4729
4730	! Read layer information (midpoint values)
4731	do iskip=1,3
4732	read (13, *)
4733	enddo
4734	do il=1,nlev-1
4735	read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
4736	n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
4737	f11(il),f12(il)
4738	pm(il)=pm(il)*100.
4739	enddo
4740	303 format (i8, 2f9.2, 10(2x,e13.7))
4741	close(13)
4742
4743	! Read aerosol layer information
4744	do iskip=1,3
4745	read (14, *)
4746	enddo
4747	read (14, '(f10.2)') aer_alpha
4748	read (14, *)
4749	read (14, *)
4750	do il=1,nlev-1
4751	read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
4752	enddo
4753	304 format (i8, f9.5, 2f8.3)
4754	close(14)
4755
4756	! Read cloud information
4757	do iskip=1,3
4758	read (15, *)
4759	enddo
4760	do il=1,nlev-1
4761	read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
4762	lwp(il)=lwp(il)/1000. ! lwp donne en g/kg
4763	iwp(il)=iwp(il)/1000. ! iwp donne en g/kg
4764	reliq(il)=reliq(il)/1000000. ! reliq donne en microns
4765	reice(il)=reice(il)/1000000. ! reice donne en microns
4766	enddo
4767	305 format (i8, f8.3, 4f9.2)
4768	close(15)
4769
4770	! Read surface albedo (weighted & unweighted) and spectral solar irradiance
4771	do iskip=1,6
4772	read (16, *)
4773	enddo
4774	do icm_1=1,ncm_1
4775	read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
4776	enddo
4777	306 format(f10.1, 2f12.5, f14.8)
4778	close(16)
4779
4780	return
4781	end subroutine read_circ
4782	!=====================================================================
4783	! Reads RTMIP input files
4784
4785	SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
4786
4787	#include "YOMCST.h"
4788
4789	real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
4790	real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
4791	integer nlev
4792
4793
4794	! Open the files
4795
4796	open (11, file='low_resolution_profile.txt', status='old')
4797
4798	! Read level information
4799	read (11, *)
4800	do il=1,nlev_rtmip
4801	read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
4802	enddo
4803	do il=1,nlev_rtmip
4804	play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb
4805	temp(il)=t(nlev_rtmip-il+1)
4806	ovap(il)=h2o(nlev_rtmip-il+1)
4807	oz(il)=o3(nlev_rtmip-il+1)
4808	enddo
4809	do il=1,39
4810	plev(il)=play(il)+(play(il+1)-play(il))/2.
4811	print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
4812	enddo
4813	plev(41)=101300.
4814	302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
4815	close(12)
4816
4817	return
4818	end subroutine read_rtmip
4819	!=====================================================================
4820
4821	! Subroutines for nudging
4822
4823	Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
4824	! ========================================================
4825	USE dimphy
4826
4827	implicit none
4828
4829	! ========================================================
4830	REAL paprs(klon,klevp1)
4831	REAL pplay(klon,klev)
4832	!
4833	! Variables d'etat
4834	REAL t(klon,klev)
4835	REAL q(klon,klev)
4836	!
4837	! Profiles cible
4838	REAL t_targ(klon,klev)
4839	REAL rh_targ(klon,klev)
4840	!
4841	INTEGER k,i
4842	REAL zx_qs
4843
4844	! Declaration des constantes et des fonctions thermodynamiques
4845	!
4846	include "YOMCST.h"
4847	include "YOETHF.h"
4848	!
4849	! ----------------------------------------
4850	! Statement functions
4851	include "FCTTRE.h"
4852	! ----------------------------------------
4853	!
4854	DO k = 1,klev
4855	DO i = 1,klon
4856	t_targ(i,k) = t(i,k)
4857	IF (t(i,k).LT.RTT) THEN
4858	zx_qs = qsats(t(i,k))/(pplay(i,k))
4859	ELSE
4860	zx_qs = qsatl(t(i,k))/(pplay(i,k))
4861	ENDIF
4862	rh_targ(i,k) = q(i,k)/zx_qs
4863	ENDDO
4864	ENDDO
4865	print *, 't_targ',t_targ
4866	print *, 'rh_targ',rh_targ
4867	!
4868	!
4869	RETURN
4870	END
4871
4872	Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
4873	! ========================================================
4874	USE dimphy
4875
4876	implicit none
4877
4878	! ========================================================
4879	REAL paprs(klon,klevp1)
4880	REAL pplay(klon,klev)
4881	!
4882	! Variables d'etat
4883	REAL u(klon,klev)
4884	REAL v(klon,klev)
4885	!
4886	! Profiles cible
4887	REAL u_targ(klon,klev)
4888	REAL v_targ(klon,klev)
4889	!
4890	INTEGER k,i
4891	!
4892	DO k = 1,klev
4893	DO i = 1,klon
4894	u_targ(i,k) = u(i,k)
4895	v_targ(i,k) = v(i,k)
4896	ENDDO
4897	ENDDO
4898	print *, 'u_targ',u_targ
4899	print *, 'v_targ',v_targ
4900	!
4901	!
4902	RETURN
4903	END
4904
4905	Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q, &
4906	& d_t,d_q)
4907	! ========================================================
4908	USE dimphy
4909
4910	implicit none
4911
4912	! ========================================================
4913	REAL dtime
4914	REAL paprs(klon,klevp1)
4915	REAL pplay(klon,klev)
4916	!
4917	! Variables d'etat
4918	REAL t(klon,klev)
4919	REAL q(klon,klev)
4920	!
4921	! Tendances
4922	REAL d_t(klon,klev)
4923	REAL d_q(klon,klev)
4924	!
4925	! Profiles cible
4926	REAL t_targ(klon,klev)
4927	REAL rh_targ(klon,klev)
4928	!
4929	! Temps de relaxation
4930	REAL tau
4931	!c DATA tau /3600./
4932	!! DATA tau /5400./
4933	DATA tau /1800./
4934	!
4935	INTEGER k,i
4936	REAL zx_qs, rh, tnew, d_rh, rhnew
4937
4938	! Declaration des constantes et des fonctions thermodynamiques
4939	!
4940	include "YOMCST.h"
4941	include "YOETHF.h"
4942	!
4943	! ----------------------------------------
4944	! Statement functions
4945	include "FCTTRE.h"
4946	! ----------------------------------------
4947	!
4948	print *,'dtime, tau ',dtime,tau
4949	print *, 't_targ',t_targ
4950	print *, 'rh_targ',rh_targ
4951	print *,'temp ',t
4952	print *,'hum ',q
4953	!
4954	DO k = 1,klev
4955	DO i = 1,klon
4956	IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
4957	IF (t(i,k).LT.RTT) THEN
4958	zx_qs = qsats(t(i,k))/(pplay(i,k))
4959	ELSE
4960	zx_qs = qsatl(t(i,k))/(pplay(i,k))
4961	ENDIF
4962	rh = q(i,k)/zx_qs
4963	!
4964	d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
4965	d_rh = 1./tau*(rh_targ(i,k)-rh)
4966	!
4967	tnew = t(i,k)+d_t(i,k)*dtime
4968	!jyg<
4969	! Formule pour q :
4970	! d_q = (1/tau) [rh_targ*qsat(T_new) - q]
4971	!
4972	! Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
4973	! qui n'etait pas correcte.
4974	!
4975	IF (tnew.LT.RTT) THEN
4976	zx_qs = qsats(tnew)/(pplay(i,k))
4977	ELSE
4978	zx_qs = qsatl(tnew)/(pplay(i,k))
4979	ENDIF
4980	!! d_q(i,k) = d_q(i,k) + d_rh*zx_qs
4981	d_q(i,k) = d_q(i,k) + (1./tau)(rh_targ(i,k)zx_qs - q(i,k))
4982	rhnew = (q(i,k)+d_q(i,k)*dtime)/zx_qs
4983	!
4984	print *,' k,d_t,rh,d_rh,rhnew,d_q ', &
4985	k,d_t(i,k),rh,d_rh,rhnew,d_q(i,k)
4986	ENDIF
4987	!
4988	ENDDO
4989	ENDDO
4990	!
4991	RETURN
4992	END
4993
4994	Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v, &
4995	& d_u,d_v)
4996	! ========================================================
4997	USE dimphy
4998
4999	implicit none
5000
5001	! ========================================================
5002	REAL dtime
5003	REAL paprs(klon,klevp1)
5004	REAL pplay(klon,klev)
5005	!
5006	! Variables d'etat
5007	REAL u(klon,klev)
5008	REAL v(klon,klev)
5009	!
5010	! Tendances
5011	REAL d_u(klon,klev)
5012	REAL d_v(klon,klev)
5013	!
5014	! Profiles cible
5015	REAL u_targ(klon,klev)
5016	REAL v_targ(klon,klev)
5017	!
5018	! Temps de relaxation
5019	REAL tau
5020	!c DATA tau /3600./
5021	DATA tau /5400./
5022	!
5023	INTEGER k,i
5024
5025	!
5026	print *,'dtime, tau ',dtime,tau
5027	print *, 'u_targ',u_targ
5028	print *, 'v_targ',v_targ
5029	print *,'zonal velocity ',u
5030	print *,'meridional velocity ',v
5031	DO k = 1,klev
5032	DO i = 1,klon
5033	IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
5034	!
5035	d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
5036	d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
5037	!
5038	print *,' k,u,d_u,v,d_v ', &
5039	k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
5040	ENDIF
5041	!
5042	ENDDO
5043	ENDDO
5044	!
5045	RETURN
5046	END
5047
5048	!=====================================================================
5049	SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas &
5050	& ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas &
5051	& ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas &
5052	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
5053	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
5054	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
5055	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas &
5056	!
5057	& ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas &
5058	& ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas &
5059	& ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas &
5060	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
5061	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas &
5062	& ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
5063
5064	implicit none
5065
5066	#include "dimensions.h"
5067
5068	!-------------------------------------------------------------------------
5069	! Vertical interpolation of generic case forcing data onto mod_casel levels
5070	!-------------------------------------------------------------------------
5071
5072	integer nlevmax
5073	parameter (nlevmax=41)
5074	integer nlev_cas,mxcalc
5075	! real play(llm), plev_prof(nlevmax)
5076	! real t_prof(nlevmax),q_prof(nlevmax)
5077	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
5078	! real ht_prof(nlevmax),vt_prof(nlevmax)
5079	! real hq_prof(nlevmax),vq_prof(nlevmax)
5080
5081	real play(llm), plev_prof_cas(nlev_cas)
5082	real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
5083	real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
5084	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
5085	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
5086	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
5087	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
5088	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
5089	real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
5090	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
5091
5092	real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
5093	real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
5094	real u_mod_cas(llm),v_mod_cas(llm)
5095	real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
5096	real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
5097	real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
5098	real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
5099	real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
5100	real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
5101
5102	integer l,k,k1,k2
5103	real frac,frac1,frac2,fact
5104
5105	do l = 1, llm
5106	print *,'debut interp2, play=',l,play(l)
5107	enddo
5108	! do l = 1, nlev_cas
5109	! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
5110	! enddo
5111
5112	do l = 1, llm
5113
5114	if (play(l).ge.plev_prof_cas(nlev_cas)) then
5115
5116	mxcalc=l
5117	print *,'debut interp2, mxcalc=',mxcalc
5118	k1=0
5119	k2=0
5120
5121	if (play(l).le.plev_prof_cas(1)) then
5122
5123	do k = 1, nlev_cas-1
5124	if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
5125	k1=k
5126	k2=k+1
5127	endif
5128	enddo
5129
5130	if (k1.eq.0 .or. k2.eq.0) then
5131	write(,) 'PB! k1, k2 = ',k1,k2
5132	write(,) 'l,play(l) = ',l,play(l)/100
5133	do k = 1, nlev_cas-1
5134	write(,) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
5135	enddo
5136	endif
5137
5138	frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
5139	t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
5140	theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
5141	thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
5142	thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
5143	qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
5144	ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
5145	qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
5146	u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
5147	v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
5148	ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
5149	vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
5150	w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
5151	omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
5152	du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
5153	hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
5154	vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
5155	dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
5156	hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
5157	vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
5158	dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
5159	ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
5160	vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
5161	dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
5162	hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
5163	vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
5164	dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
5165	hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
5166	vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
5167	dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
5168
5169	else !play>plev_prof_cas(1)
5170
5171	k1=1
5172	k2=2
5173	print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
5174	frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
5175	frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
5176	t_mod_cas(l)= frac1t_prof_cas(k1) - frac2t_prof_cas(k2)
5177	theta_mod_cas(l)= frac1th_prof_cas(k1) - frac2th_prof_cas(k2)
5178	thv_mod_cas(l)= frac1thv_prof_cas(k1) - frac2thv_prof_cas(k2)
5179	thl_mod_cas(l)= frac1thl_prof_cas(k1) - frac2thl_prof_cas(k2)
5180	qv_mod_cas(l)= frac1qv_prof_cas(k1) - frac2qv_prof_cas(k2)
5181	ql_mod_cas(l)= frac1ql_prof_cas(k1) - frac2ql_prof_cas(k2)
5182	qi_mod_cas(l)= frac1qi_prof_cas(k1) - frac2qi_prof_cas(k2)
5183	u_mod_cas(l)= frac1u_prof_cas(k1) - frac2u_prof_cas(k2)
5184	v_mod_cas(l)= frac1v_prof_cas(k1) - frac2v_prof_cas(k2)
5185	ug_mod_cas(l)= frac1ug_prof_cas(k1) - frac2ug_prof_cas(k2)
5186	vg_mod_cas(l)= frac1vg_prof_cas(k1) - frac2vg_prof_cas(k2)
5187	w_mod_cas(l)= frac1vitw_prof_cas(k1) - frac2vitw_prof_cas(k2)
5188	omega_mod_cas(l)= frac1omega_prof_cas(k1) - frac2omega_prof_cas(k2)
5189	du_mod_cas(l)= frac1du_prof_cas(k1) - frac2du_prof_cas(k2)
5190	hu_mod_cas(l)= frac1hu_prof_cas(k1) - frac2hu_prof_cas(k2)
5191	vu_mod_cas(l)= frac1vu_prof_cas(k1) - frac2vu_prof_cas(k2)
5192	dv_mod_cas(l)= frac1dv_prof_cas(k1) - frac2dv_prof_cas(k2)
5193	hv_mod_cas(l)= frac1hv_prof_cas(k1) - frac2hv_prof_cas(k2)
5194	vv_mod_cas(l)= frac1vv_prof_cas(k1) - frac2vv_prof_cas(k2)
5195	dt_mod_cas(l)= frac1dt_prof_cas(k1) - frac2dt_prof_cas(k2)
5196	ht_mod_cas(l)= frac1ht_prof_cas(k1) - frac2ht_prof_cas(k2)
5197	vt_mod_cas(l)= frac1vt_prof_cas(k1) - frac2vt_prof_cas(k2)
5198	dth_mod_cas(l)= frac1dth_prof_cas(k1) - frac2dth_prof_cas(k2)
5199	hth_mod_cas(l)= frac1hth_prof_cas(k1) - frac2hth_prof_cas(k2)
5200	vth_mod_cas(l)= frac1vth_prof_cas(k1) - frac2vth_prof_cas(k2)
5201	dq_mod_cas(l)= frac1dq_prof_cas(k1) - frac2dq_prof_cas(k2)
5202	hq_mod_cas(l)= frac1hq_prof_cas(k1) - frac2hq_prof_cas(k2)
5203	vq_mod_cas(l)= frac1vq_prof_cas(k1) - frac2vq_prof_cas(k2)
5204	dtrad_mod_cas(l)= frac1dtrad_prof_cas(k1) - frac2dtrad_prof_cas(k2)
5205
5206	endif ! play.le.plev_prof_cas(1)
5207
5208	else ! above max altitude of forcing file
5209
5210	!jyg
5211	fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
5212	fact = max(fact,0.) !jyg
5213	fact = exp(-fact) !jyg
5214	t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg
5215	theta_mod_cas(l)= th_prof_cas(nlev_cas) !jyg
5216	thv_mod_cas(l)= thv_prof_cas(nlev_cas) !jyg
5217	thl_mod_cas(l)= thl_prof_cas(nlev_cas) !jyg
5218	qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact !jyg
5219	ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact !jyg
5220	qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact !jyg
5221	u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg
5222	v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg
5223	ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg
5224	vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg
5225	w_mod_cas(l)= 0.0 !jyg
5226	du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
5227	hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg
5228	vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg
5229	dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
5230	hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg
5231	vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg
5232	dt_mod_cas(l)= dt_prof_cas(nlev_cas)
5233	ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg
5234	vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg
5235	dth_mod_cas(l)= dth_prof_cas(nlev_cas)
5236	hth_mod_cas(l)= hth_prof_cas(nlev_cas) !jyg
5237	vth_mod_cas(l)= vth_prof_cas(nlev_cas) !jyg
5238	dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
5239	hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg
5240	vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg
5241	dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg
5242
5243	endif ! play
5244
5245	enddo ! l
5246
5247	! do l = 1,llm
5248	! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
5249	! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
5250	! enddo
5251
5252	return
5253	end
5254	!*****************************************************************************
5255
5256

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

Download in other formats: