Context Navigation

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

Last change on this file since 2311 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM

File size: 143.3 KB

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

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

Download in other formats: