Context Navigation

source: LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 2672

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

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

Download in other formats: