Context Navigation

source: LMDZ5/branches/IPSLCM6.0.11/libf/phylmd/dyn1d/1DUTILS.h @ 5416

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

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

Download in other formats: