Context Navigation

source: LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 3518

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

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

Download in other formats: