Context Navigation

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

Last change on this file since 3957 was 2857, checked in by jyg, 8 years ago
Bug fix concerning 1D restart
File size: 173.8 KB

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

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

Download in other formats: