Context Navigation

source: LMDZ5/branches/testing/libf/phylmd/dyn1d/1DUTILS.h @ 2471

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

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

Download in other formats: