Context Navigation

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

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1997:2055 into testing branch
File size: 111.6 KB

Line
1	#include "../dyn3d/conf_gcm.F"
2	#include "../dyn3d_common/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	IMPLICIT NONE
18	!-----------------------------------------------------------------------
19	! Auteurs : A. Lahellec .
20	!
21	! Declarations :
22	! --------------
23
24	#include "compar1d.h"
25	#include "flux_arp.h"
26	#include "tsoilnudge.h"
27	#include "fcg_gcssold.h"
28	#include "fcg_racmo.h"
29	#include "iniprint.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	! = 40 ==> forcing_GCSSold = .true.
102	! initial profile from GCSS file
103	! LS convergence imposed from GCSS file
104	! = 50 ==> forcing_fire = .true.
105	! = 59 ==> forcing_sandu = .true.
106	! initial profiles from sanduref file: see prof.inp.001
107	! SST varying with time and divergence constante: see ifa_sanduref.txt file
108	! Radiation has to be computed interactively
109	! = 60 ==> forcing_astex = .true.
110	! initial profiles from file: see prof.inp.001
111	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
112	! Radiation has to be computed interactively
113	! = 61 ==> forcing_armcu = .true.
114	! initial profiles from file: see prof.inp.001
115	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
116	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
117	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
118	! Radiation to be switched off
119	!
120	forcing_type = 0
121	CALL getin('forcing_type',forcing_type)
122	imp_fcg_gcssold = .false.
123	ts_fcg_gcssold = .false.
124	Tp_fcg_gcssold = .false.
125	Tp_ini_gcssold = .false.
126	xTurb_fcg_gcssold = .false.
127	IF (forcing_type .eq.40) THEN
128	CALL getin('imp_fcg',imp_fcg_gcssold)
129	CALL getin('ts_fcg',ts_fcg_gcssold)
130	CALL getin('tp_fcg',Tp_fcg_gcssold)
131	CALL getin('tp_ini',Tp_ini_gcssold)
132	CALL getin('turb_fcg',xTurb_fcg_gcssold)
133	ENDIF
134
135	!Config Key = ok_flux_surf
136	!Config Desc = forcage ou non par les flux de surface
137	!Config Def = false
138	!Config Help = forcage ou non par les flux de surface
139	ok_flux_surf =.false.
140	CALL getin('ok_flux_surf',ok_flux_surf)
141
142	!Config Key = ok_old_disvert
143	!Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
144	!Config Def = false
145	!Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
146	ok_old_disvert = .FALSE.
147	CALL getin('ok_old_disvert',ok_old_disvert)
148
149	!Config Key = time_ini
150	!Config Desc = meaningless in this case
151	!Config Def = 0.
152	!Config Help =
153	tsurf = 0.
154	CALL getin('time_ini',time_ini)
155
156	!Config Key = rlat et rlon
157	!Config Desc = latitude et longitude
158	!Config Def = 0.0 0.0
159	!Config Help = fixe la position de la colonne
160	xlat = 0.
161	xlon = 0.
162	CALL getin('rlat',xlat)
163	CALL getin('rlon',xlon)
164
165	!Config Key = airephy
166	!Config Desc = Grid cell area
167	!Config Def = 1.e11
168	!Config Help =
169	airefi = 1.e11
170	CALL getin('airephy',airefi)
171
172	!Config Key = nat_surf
173	!Config Desc = surface type
174	!Config Def = 0 (ocean)
175	!Config Help = 0=ocean,1=land,2=glacier,3=banquise
176	nat_surf = 0.
177	CALL getin('nat_surf',nat_surf)
178
179	!Config Key = tsurf
180	!Config Desc = surface temperature
181	!Config Def = 290.
182	!Config Help = not used if type_ts_forcing=1 in lmdz1d.F
183	tsurf = 290.
184	CALL getin('tsurf',tsurf)
185
186	!Config Key = psurf
187	!Config Desc = surface pressure
188	!Config Def = 102400.
189	!Config Help =
190	psurf = 102400.
191	CALL getin('psurf',psurf)
192
193	!Config Key = zsurf
194	!Config Desc = surface altitude
195	!Config Def = 0.
196	!Config Help =
197	zsurf = 0.
198	CALL getin('zsurf',zsurf)
199
200	!Config Key = rugos
201	!Config Desc = coefficient de frottement
202	!Config Def = 0.0001
203	!Config Help = calcul du Cdrag
204	rugos = 0.0001
205	CALL getin('rugos',rugos)
206
207	!Config Key = wtsurf et wqsurf
208	!Config Desc = ???
209	!Config Def = 0.0 0.0
210	!Config Help =
211	wtsurf = 0.0
212	wqsurf = 0.0
213	CALL getin('wtsurf',wtsurf)
214	CALL getin('wqsurf',wqsurf)
215
216	!Config Key = albedo
217	!Config Desc = albedo
218	!Config Def = 0.09
219	!Config Help =
220	albedo = 0.09
221	CALL getin('albedo',albedo)
222
223	!Config Key = agesno
224	!Config Desc = age de la neige
225	!Config Def = 30.0
226	!Config Help =
227	xagesno = 30.0
228	CALL getin('agesno',xagesno)
229
230	!Config Key = restart_runoff
231	!Config Desc = age de la neige
232	!Config Def = 30.0
233	!Config Help =
234	restart_runoff = 0.0
235	CALL getin('restart_runoff',restart_runoff)
236
237	!Config Key = qsolinp
238	!Config Desc = initial bucket water content (kg/m2) when land (5std)
239	!Config Def = 30.0
240	!Config Help =
241	qsolinp = 1.
242	CALL getin('qsolinp',qsolinp)
243
244	!Config Key = zpicinp
245	!Config Desc = denivellation orographie
246	!Config Def = 300.
247	!Config Help = input brise
248	zpicinp = 300.
249	CALL getin('zpicinp',zpicinp)
250	!Config key = nudge_tsoil
251	!Config Desc = activation of soil temperature nudging
252	!Config Def = .FALSE.
253	!Config Help = ...
254
255	nudge_tsoil=.FALSE.
256	CALL getin('nudge_tsoil',nudge_tsoil)
257
258	!Config key = isoil_nudge
259	!Config Desc = level number where soil temperature is nudged
260	!Config Def = 3
261	!Config Help = ...
262
263	isoil_nudge=3
264	CALL getin('isoil_nudge',isoil_nudge)
265
266	!Config key = Tsoil_nudge
267	!Config Desc = target temperature for tsoil(isoil_nudge)
268	!Config Def = 300.
269	!Config Help = ...
270
271	Tsoil_nudge=300.
272	CALL getin('Tsoil_nudge',Tsoil_nudge)
273
274	!Config key = tau_soil_nudge
275	!Config Desc = nudging relaxation time for tsoil
276	!Config Def = 3600.
277	!Config Help = ...
278
279	tau_soil_nudge=3600.
280	CALL getin('tau_soil_nudge',tau_soil_nudge)
281
282
283
284
285	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
286	write(lunout,*)' Configuration des parametres du gcm1D: '
287	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
288	write(lunout,*)' restart = ', restart
289	write(lunout,*)' forcing_type = ', forcing_type
290	write(lunout,*)' time_ini = ', time_ini
291	write(lunout,*)' rlat = ', xlat
292	write(lunout,*)' rlon = ', xlon
293	write(lunout,*)' airephy = ', airefi
294	write(lunout,*)' nat_surf = ', nat_surf
295	write(lunout,*)' tsurf = ', tsurf
296	write(lunout,*)' psurf = ', psurf
297	write(lunout,*)' zsurf = ', zsurf
298	write(lunout,*)' rugos = ', rugos
299	write(lunout,*)' wtsurf = ', wtsurf
300	write(lunout,*)' wqsurf = ', wqsurf
301	write(lunout,*)' albedo = ', albedo
302	write(lunout,*)' xagesno = ', xagesno
303	write(lunout,*)' restart_runoff = ', restart_runoff
304	write(lunout,*)' qsolinp = ', qsolinp
305	write(lunout,*)' zpicinp = ', zpicinp
306	write(lunout,*)' nudge_tsoil = ', nudge_tsoil
307	write(lunout,*)' isoil_nudge = ', isoil_nudge
308	write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
309	write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
310	IF (forcing_type .eq.40) THEN
311	write(lunout,*) '--- Forcing type GCSS Old --- with:'
312	write(lunout,*)'imp_fcg',imp_fcg_gcssold
313	write(lunout,*)'ts_fcg',ts_fcg_gcssold
314	write(lunout,*)'tp_fcg',Tp_fcg_gcssold
315	write(lunout,*)'tp_ini',Tp_ini_gcssold
316	write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
317	ENDIF
318
319	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
320	write(lunout,*)
321	!
322	RETURN
323	END
324	!
325	! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
326	!
327	!
328	SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs, &
329	& ucov,vcov,temp,q,omega2)
330	USE dimphy
331	USE mod_grid_phy_lmdz
332	USE mod_phys_lmdz_para
333	USE iophy
334	USE phys_state_var_mod
335	USE iostart
336	USE write_field_phy
337	USE infotrac
338	use control_mod
339
340	IMPLICIT NONE
341	!=======================================================
342	! Ecriture du fichier de redemarrage sous format NetCDF
343	!=======================================================
344	! Declarations:
345	! -------------
346	#include "dimensions.h"
347	#include "comconst.h"
348	#include "temps.h"
349	!!#include "control.h"
350	#include "logic.h"
351	#include "netcdf.inc"
352
353	! Arguments:
354	! ----------
355	CHARACTER() fichnom
356	!Al1 plev tronque pour .nc mais plev(klev+1):=0
357	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
358	real :: presnivs(klon,klev)
359	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
360	real :: q(klon,klev,nqtot),omega2(klon,klev)
361	! real :: ug(klev),vg(klev),fcoriolis
362	real :: phis(klon)
363
364	! Variables locales pour NetCDF:
365	! ------------------------------
366	INTEGER iq
367	INTEGER length
368	PARAMETER (length = 100)
369	REAL tab_cntrl(length) ! tableau des parametres du run
370	character*4 nmq(nqtot)
371	character*12 modname
372	character*80 abort_message
373	LOGICAL found
374
375	modname = 'dyn1deta0 : '
376	nmq(1)="vap"
377	nmq(2)="cond"
378	do iq=3,nqtot
379	write(nmq(iq),'("tra",i1)') iq-2
380	enddo
381	print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
382	CALL open_startphy(fichnom)
383	print*,'after open startphy ',fichnom,nmq
384
385	!
386	! Lecture des parametres de controle:
387	!
388	CALL get_var("controle",tab_cntrl)
389
390
391	im = tab_cntrl(1)
392	jm = tab_cntrl(2)
393	lllm = tab_cntrl(3)
394	day_ref = tab_cntrl(4)
395	annee_ref = tab_cntrl(5)
396	! rad = tab_cntrl(6)
397	! omeg = tab_cntrl(7)
398	! g = tab_cntrl(8)
399	! cpp = tab_cntrl(9)
400	! kappa = tab_cntrl(10)
401	! daysec = tab_cntrl(11)
402	! dtvr = tab_cntrl(12)
403	! etot0 = tab_cntrl(13)
404	! ptot0 = tab_cntrl(14)
405	! ztot0 = tab_cntrl(15)
406	! stot0 = tab_cntrl(16)
407	! ang0 = tab_cntrl(17)
408	! pa = tab_cntrl(18)
409	! preff = tab_cntrl(19)
410	!
411	! clon = tab_cntrl(20)
412	! clat = tab_cntrl(21)
413	! grossismx = tab_cntrl(22)
414	! grossismy = tab_cntrl(23)
415	!
416	IF ( tab_cntrl(24).EQ.1. ) THEN
417	fxyhypb =.true.
418	! dzoomx = tab_cntrl(25)
419	! dzoomy = tab_cntrl(26)
420	! taux = tab_cntrl(28)
421	! tauy = tab_cntrl(29)
422	ELSE
423	fxyhypb = .false.
424	ysinus = .false.
425	IF( tab_cntrl(27).EQ.1. ) ysinus =.true.
426	ENDIF
427
428	day_ini = tab_cntrl(30)
429	itau_dyn = tab_cntrl(31)
430	! .................................................................
431	!
432	!
433	! PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
434	!Al1
435	Print*,'day_ref,annee_ref,day_ini,itau_dyn', &
436	& day_ref,annee_ref,day_ini,itau_dyn
437
438	! Lecture des champs
439	!
440	plev(1,klev+1)=0.
441	CALL get_field("plev",plev,found)
442	IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
443	CALL get_field("play",play,found)
444	IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
445	CALL get_field("phi",phi,found)
446	IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
447	CALL get_field("phis",phis,found)
448	IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
449	CALL get_field("presnivs",presnivs,found)
450	IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
451	CALL get_field("ucov",ucov,found)
452	IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
453	CALL get_field("vcov",vcov,found)
454	IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
455	CALL get_field("temp",temp,found)
456	IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
457	CALL get_field("omega2",omega2,found)
458	IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
459
460	Do iq=1,nqtot
461	CALL get_field("q"//nmq(iq),q(:,:,iq),found)
462	IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
463	EndDo
464
465	CALL close_startphy
466	print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
467	!
468	RETURN
469	END
470	!
471	! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
472	!
473	!
474	SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs, &
475	& ucov,vcov,temp,q,omega2)
476	USE dimphy
477	USE mod_grid_phy_lmdz
478	USE mod_phys_lmdz_para
479	USE phys_state_var_mod
480	USE iostart
481	USE infotrac
482	use control_mod
483
484	IMPLICIT NONE
485	!=======================================================
486	! Ecriture du fichier de redemarrage sous format NetCDF
487	!=======================================================
488	! Declarations:
489	! -------------
490	#include "dimensions.h"
491	#include "comconst.h"
492	#include "temps.h"
493	!!#include "control.h"
494	#include "logic.h"
495	#include "netcdf.inc"
496
497	! Arguments:
498	! ----------
499	CHARACTER() fichnom
500	!Al1 plev tronque pour .nc mais plev(klev+1):=0
501	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
502	real :: presnivs(klon,klev)
503	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
504	real :: q(klon,klev,nqtot)
505	real :: omega2(klon,klev),rho(klon,klev+1)
506	! real :: ug(klev),vg(klev),fcoriolis
507	real :: phis(klon)
508
509	! Variables locales pour NetCDF:
510	! ------------------------------
511	INTEGER nid
512	INTEGER ierr
513	INTEGER iq,l
514	INTEGER length
515	PARAMETER (length = 100)
516	REAL tab_cntrl(length) ! tableau des parametres du run
517	character*4 nmq(nqtot)
518	character*20 modname
519	character*80 abort_message
520	!
521	INTEGER nb
522	SAVE nb
523	DATA nb / 0 /
524
525	CALL open_restartphy(fichnom)
526	print*,'redm1 ',fichnom,klon,klev,nqtot
527	nmq(1)="vap"
528	nmq(2)="cond"
529	nmq(3)="tra1"
530	nmq(4)="tra2"
531
532	modname = 'dyn1dredem'
533	ierr = NF_OPEN(fichnom, NF_WRITE, nid)
534	IF (ierr .NE. NF_NOERR) THEN
535	abort_message="Pb. d ouverture "//fichnom
536	CALL abort_gcm('Modele 1D',abort_message,1)
537	ENDIF
538
539	DO l=1,length
540	tab_cntrl(l) = 0.
541	ENDDO
542	tab_cntrl(1) = FLOAT(iim)
543	tab_cntrl(2) = FLOAT(jjm)
544	tab_cntrl(3) = FLOAT(llm)
545	tab_cntrl(4) = FLOAT(day_ref)
546	tab_cntrl(5) = FLOAT(annee_ref)
547	tab_cntrl(6) = rad
548	tab_cntrl(7) = omeg
549	tab_cntrl(8) = g
550	tab_cntrl(9) = cpp
551	tab_cntrl(10) = kappa
552	tab_cntrl(11) = daysec
553	tab_cntrl(12) = dtvr
554	! tab_cntrl(13) = etot0
555	! tab_cntrl(14) = ptot0
556	! tab_cntrl(15) = ztot0
557	! tab_cntrl(16) = stot0
558	! tab_cntrl(17) = ang0
559	! tab_cntrl(18) = pa
560	! tab_cntrl(19) = preff
561	!
562	! ..... parametres pour le zoom ......
563
564	! tab_cntrl(20) = clon
565	! tab_cntrl(21) = clat
566	! tab_cntrl(22) = grossismx
567	! tab_cntrl(23) = grossismy
568	!
569	IF ( fxyhypb ) THEN
570	tab_cntrl(24) = 1.
571	! tab_cntrl(25) = dzoomx
572	! tab_cntrl(26) = dzoomy
573	tab_cntrl(27) = 0.
574	! tab_cntrl(28) = taux
575	! tab_cntrl(29) = tauy
576	ELSE
577	tab_cntrl(24) = 0.
578	! tab_cntrl(25) = dzoomx
579	! tab_cntrl(26) = dzoomy
580	tab_cntrl(27) = 0.
581	tab_cntrl(28) = 0.
582	tab_cntrl(29) = 0.
583	IF( ysinus ) tab_cntrl(27) = 1.
584	ENDIF
585	!Al1 iday_end -> day_end
586	tab_cntrl(30) = FLOAT(day_end)
587	tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
588	!
589	CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
590	!
591
592	! Ecriture/extension de la coordonnee temps
593
594	nb = nb + 1
595
596	! Ecriture des champs
597	!
598	CALL put_field("plev","p interfaces sauf la nulle",plev)
599	CALL put_field("play","",play)
600	CALL put_field("phi","geopotentielle",phi)
601	CALL put_field("phis","geopotentiell de surface",phis)
602	CALL put_field("presnivs","",presnivs)
603	CALL put_field("ucov","",ucov)
604	CALL put_field("vcov","",vcov)
605	CALL put_field("temp","",temp)
606	CALL put_field("omega2","",omega2)
607
608	Do iq=1,nqtot
609	CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs", &
610	& q(:,:,iq))
611	EndDo
612	CALL close_restartphy
613
614	!
615	RETURN
616	END
617	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
618	IMPLICIT NONE
619	!=======================================================================
620	! passage d'un champ de la grille scalaire a la grille physique
621	!=======================================================================
622
623	!-----------------------------------------------------------------------
624	! declarations:
625	! -------------
626
627	INTEGER im,jm,ngrid,nfield
628	REAL pdyn(im,jm,nfield)
629	REAL pfi(ngrid,nfield)
630
631	INTEGER i,j,ifield,ig
632
633	!-----------------------------------------------------------------------
634	! calcul:
635	! -------
636
637	DO ifield=1,nfield
638	! traitement des poles
639	DO i=1,im
640	pdyn(i,1,ifield)=pfi(1,ifield)
641	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
642	ENDDO
643
644	! traitement des point normaux
645	DO j=2,jm-1
646	ig=2+(j-2)*(im-1)
647	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
648	pdyn(im,j,ifield)=pdyn(1,j,ifield)
649	ENDDO
650	ENDDO
651
652	RETURN
653	END
654
655
656
657	SUBROUTINE abort_gcm(modname, message, ierr)
658
659	USE IOIPSL
660	!
661	! Stops the simulation cleanly, closing files and printing various
662	! comments
663	!
664	! Input: modname = name of calling program
665	! message = stuff to print
666	! ierr = severity of situation ( = 0 normal )
667
668	character(len=*) modname
669	integer ierr
670	character(len=*) message
671
672	write(,) 'in abort_gcm'
673	call histclo
674	! call histclo(2)
675	! call histclo(3)
676	! call histclo(4)
677	! call histclo(5)
678	write(,) 'out of histclo'
679	write(,) 'Stopping in ', modname
680	write(,) 'Reason = ',message
681	call getin_dump
682	!
683	if (ierr .eq. 0) then
684	write(,) 'Everything is cool'
685	else
686	write(,) 'Houston, we have a problem ', ierr
687	endif
688	STOP
689	END
690	REAL FUNCTION fq_sat(kelvin, millibar)
691	!
692	IMPLICIT none
693	!======================================================================
694	! Autheur(s): Z.X. Li (LMD/CNRS)
695	! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
696	!======================================================================
697	! Arguments:
698	! kelvin---input-R: temperature en Kelvin
699	! millibar--input-R: pression en mb
700	!
701	! fq_sat----output-R: vapeur d'eau saturante en kg/kg
702	!======================================================================
703	!
704	REAL kelvin, millibar
705	!
706	REAL r2es
707	PARAMETER (r2es=611.14 *18.0153/28.9644)
708	!
709	REAL r3les, r3ies, r3es
710	PARAMETER (R3LES=17.269)
711	PARAMETER (R3IES=21.875)
712	!
713	REAL r4les, r4ies, r4es
714	PARAMETER (R4LES=35.86)
715	PARAMETER (R4IES=7.66)
716	!
717	REAL rtt
718	PARAMETER (rtt=273.16)
719	!
720	REAL retv
721	PARAMETER (retv=28.9644/18.0153 - 1.0)
722	!
723	REAL zqsat
724	REAL temp, pres
725	! ------------------------------------------------------------------
726	!
727	!
728	temp = kelvin
729	pres = millibar * 100.0
730	! write(,)'kelvin,millibar=',kelvin,millibar
731	! write(,)'temp,pres=',temp,pres
732	!
733	IF (temp .LE. rtt) THEN
734	r3es = r3ies
735	r4es = r4ies
736	ELSE
737	r3es = r3les
738	r4es = r4les
739	ENDIF
740	!
741	zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
742	zqsat=MIN(0.5,ZQSAT)
743	zqsat=zqsat/(1.-retv *zqsat)
744	!
745	fq_sat = zqsat
746	!
747	RETURN
748	END
749
750	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
751	IMPLICIT NONE
752	!=======================================================================
753	! passage d'un champ de la grille scalaire a la grille physique
754	!=======================================================================
755
756	!-----------------------------------------------------------------------
757	! declarations:
758	! -------------
759
760	INTEGER im,jm,ngrid,nfield
761	REAL pdyn(im,jm,nfield)
762	REAL pfi(ngrid,nfield)
763
764	INTEGER j,ifield,ig
765
766	!-----------------------------------------------------------------------
767	! calcul:
768	! -------
769
770	IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1) &
771	& STOP 'probleme de dim'
772	! traitement des poles
773	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
774	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
775
776	! traitement des point normaux
777	DO ifield=1,nfield
778	DO j=2,jm-1
779	ig=2+(j-2)*(im-1)
780	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
781	ENDDO
782	ENDDO
783
784	RETURN
785	END
786
787	SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
788
789	! Ancienne version disvert dont on a modifie nom pour utiliser
790	! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
791	! (MPL 18092012)
792	!
793	! Auteur : P. Le Van .
794	!
795	IMPLICIT NONE
796
797	#include "dimensions.h"
798	#include "paramet.h"
799	!
800	!=======================================================================
801	!
802	!
803	! s = sigma ** kappa : coordonnee verticale
804	! dsig(l) : epaisseur de la couche l ds la coord. s
805	! sig(l) : sigma a l'interface des couches l et l-1
806	! ds(l) : distance entre les couches l et l-1 en coord.s
807	!
808	!=======================================================================
809	!
810	REAL pa,preff
811	REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
812	REAL presnivs(llm)
813	!
814	! declarations:
815	! -------------
816	!
817	REAL sig(llm+1),dsig(llm)
818	!
819	INTEGER l
820	REAL snorm
821	REAL alpha,beta,gama,delta,deltaz,h
822	INTEGER np,ierr
823	REAL pi,x
824	REAL :: scaleheight=8.
825
826	!-----------------------------------------------------------------------
827	!
828	pi=2.*ASIN(1.)
829
830	OPEN(99,file='sigma.def',status='old',form='formatted', &
831	& iostat=ierr)
832
833	!-----------------------------------------------------------------------
834	! cas 1 on lit les options dans sigma.def:
835	! ----------------------------------------
836
837	IF (ierr.eq.0) THEN
838
839	print*,'WARNING!!! on lit les options dans sigma.def'
840	READ(99,*) deltaz
841	READ(99,*) h
842	READ(99,*) beta
843	READ(99,*) gama
844	READ(99,*) delta
845	READ(99,*) np
846	CLOSE(99)
847	alpha=deltaz/(llm*h)
848	!
849
850	DO 1 l = 1, llm
851	dsig(l) = (alpha+(1.-alpha)exp(-beta(llm-l)))* &
852	& ( (tanh(gamal)/tanh(gamallm))**np + &
853	& (1.-l/FLOAT(llm))*delta )
854	1 CONTINUE
855
856	sig(1)=1.
857	DO 101 l=1,llm-1
858	sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
859	101 CONTINUE
860	sig(llm+1)=0.
861
862	DO 2 l = 1, llm
863	dsig(l) = sig(l)-sig(l+1)
864	2 CONTINUE
865	!
866
867	ELSE
868	!-----------------------------------------------------------------------
869	! cas 2 ancienne discretisation (LMD5...):
870	! ----------------------------------------
871
872	PRINT*,'WARNING!!! Ancienne discretisation verticale'
873
874	h=7.
875	snorm = 0.
876	DO l = 1, llm
877	x = 2.asin(1.) (FLOAT(l)-0.5) / float(llm+1)
878	dsig(l) = 1.0 + 7.0 * SIN(x)**2
879	snorm = snorm + dsig(l)
880	ENDDO
881	snorm = 1./snorm
882	DO l = 1, llm
883	dsig(l) = dsig(l)*snorm
884	ENDDO
885	sig(llm+1) = 0.
886	DO l = llm, 1, -1
887	sig(l) = sig(l+1) + dsig(l)
888	ENDDO
889
890	ENDIF
891
892
893	DO l=1,llm
894	nivsigs(l) = FLOAT(l)
895	ENDDO
896
897	DO l=1,llmp1
898	nivsig(l)= FLOAT(l)
899	ENDDO
900
901	!
902	! .... Calculs de ap(l) et de bp(l) ....
903	! .........................................
904	!
905	!
906	! ..... pa et preff sont lus sur les fichiers start par lectba .....
907	!
908
909	bp(llmp1) = 0.
910
911	DO l = 1, llm
912	!c
913	!cc ap(l) = 0.
914	!cc bp(l) = sig(l)
915
916	bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
917	ap(l) = pa * ( sig(l) - bp(l) )
918	!
919	ENDDO
920	ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
921
922	PRINT *,' BP '
923	PRINT *, bp
924	PRINT *,' AP '
925	PRINT *, ap
926
927	print*,'scaleheight=',scaleheight
928	DO l = 1, llm
929	dpres(l) = bp(l) - bp(l+1)
930	presnivs(l) = 0.5 ( ap(l)+bp(l)preff + ap(l+1)+bp(l+1)*preff )
931	write(, )'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
932	log(preff/presnivs(l))*scaleheight &
933	, ' DZ ~ ', scaleheightlog((ap(l)+bp(l)preff)/ &
934	max(ap(l+1)+bp(l+1)*preff, 1.e-10))
935	ENDDO
936
937
938	PRINT *,' PRESNIVS '
939	PRINT *,presnivs
940
941	RETURN
942	END
943
944	!======================================================================
945	SUBROUTINE read_tsurf1d(knon,sst_out)
946
947	! This subroutine specifies the surface temperature to be used in 1D simulations
948
949	USE dimphy, ONLY : klon
950
951	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
952	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
953
954	INTEGER :: i
955	! COMMON defined in lmdz1d.F:
956	real ts_cur
957	common /sst_forcing/ts_cur
958
959	DO i = 1, knon
960	sst_out(i) = ts_cur
961	ENDDO
962
963	END SUBROUTINE read_tsurf1d
964
965	!===============================================================
966	subroutine advect_vert(llm,w,dt,q,plev)
967	!===============================================================
968	! Schema amont pour l'advection verticale en 1D
969	! w est la vitesse verticale dp/dt en Pa/s
970	! Traitement en volumes finis
971	! d / dt ( zm q ) = delta_z ( omega q )
972	! d / dt ( zm ) = delta_z ( omega )
973	! avec zm = delta_z ( p )
974	! si * designe la valeur au pas de temps t+dt
975	! zm(l) q(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
976	! zm*(l) -zm(l) = w(l+1) - w(l)
977	! avec w=omega * dt
978	!---------------------------------------------------------------
979	implicit none
980	! arguments
981	integer llm
982	real w(llm+1),q(llm),plev(llm+1),dt
983
984	! local
985	integer l
986	real zwq(llm+1),zm(llm+1),zw(llm+1)
987	real qold
988
989	!---------------------------------------------------------------
990
991	do l=1,llm
992	zw(l)=dt*w(l)
993	zm(l)=plev(l)-plev(l+1)
994	zwq(l)=q(l)*zw(l)
995	enddo
996	zwq(llm+1)=0.
997	zw(llm+1)=0.
998
999	do l=1,llm
1000	qold=q(l)
1001	q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1002	print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1003	enddo
1004
1005
1006	return
1007	end
1008
1009	!===============================================================
1010
1011
1012	SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, &
1013	& q,temp,u,v,play)
1014	!itlmd
1015	!----------------------------------------------------------------------
1016	! Calcul de l'advection verticale (ascendance et subsidence) de
1017	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1018	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1019	! sans WTG rajouter une advection horizontale
1020	!----------------------------------------------------------------------
1021	implicit none
1022	#include "YOMCST.h"
1023	! argument
1024	integer llm
1025	real omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1026	real d_u_va(llm), d_v_va(llm)
1027	real q(llm,3),temp(llm)
1028	real u(llm),v(llm)
1029	real play(llm)
1030	! interne
1031	integer l
1032	real alpha,omgdown,omgup
1033
1034	do l= 1,llm
1035	if(l.eq.1) then
1036	!si omgup pour la couche 1, alors tendance nulle
1037	omgdown=max(omega(2),0.0)
1038	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1039	d_t_va(l)= alpha(omgdown)-omgdown(temp(l)-temp(l+1)) &
1040	& /(play(l)-play(l+1))
1041
1042	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))
1043
1044	d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))
1045	d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))
1046
1047
1048	elseif(l.eq.llm) then
1049	omgup=min(omega(l),0.0)
1050	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1051	d_t_va(l)= alpha*(omgup)- &
1052
1053	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1054	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1055	d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1056	d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1057	d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1058
1059	else
1060	omgup=min(omega(l),0.0)
1061	omgdown=max(omega(l+1),0.0)
1062	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
1063	d_t_va(l)= alpha(omgup+omgdown)-omgdown(temp(l)-temp(l+1)) &
1064	& /(play(l)-play(l+1))- &
1065	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1066	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1067	! print*, ' ??? '
1068
1069	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:)) &
1070	& /(play(l)-play(l+1))- &
1071	& omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1072	d_u_va(l)= -omgdown*(u(l)-u(l+1)) &
1073	& /(play(l)-play(l+1))- &
1074	& omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1075	d_v_va(l)= -omgdown*(v(l)-v(l+1)) &
1076	& /(play(l)-play(l+1))- &
1077	& omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1078
1079	endif
1080
1081	enddo
1082	!fin itlmd
1083	return
1084	end
1085	! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1086	SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va, &
1087	& q,temp,u,v,play)
1088	!itlmd
1089	!----------------------------------------------------------------------
1090	! Calcul de l'advection verticale (ascendance et subsidence) de
1091	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1092	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1093	! sans WTG rajouter une advection horizontale
1094	!----------------------------------------------------------------------
1095	implicit none
1096	#include "YOMCST.h"
1097	! argument
1098	integer llm,nqtot
1099	real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1100	! real d_u_va(llm), d_v_va(llm)
1101	real q(llm,nqtot),temp(llm)
1102	real u(llm),v(llm)
1103	real play(llm)
1104	real cor(llm)
1105	! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1106	real dph(llm),dqdp(llm),dtdp(llm)
1107	! interne
1108	integer k
1109	real omdn,omup
1110
1111	! dudp=0.
1112	! dvdp=0.
1113	dqdp=0.
1114	dtdp=0.
1115	! d_u_va=0.
1116	! d_v_va=0.
1117
1118	cor(:) = rkappatemp(1.+q(:,1)rv/rd)/(play(1.+q(:,1)))
1119
1120
1121	do k=2,llm-1
1122
1123	dph (k-1) = (play(k )- play(k-1 ))
1124	! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
1125	! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1)
1126	dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1)
1127	dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1)
1128
1129	enddo
1130
1131	! dudp ( llm ) = dudp ( llm-1 )
1132	! dvdp ( llm ) = dvdp ( llm-1 )
1133	dqdp ( llm ) = dqdp ( llm-1 )
1134	dtdp ( llm ) = dtdp ( llm-1 )
1135
1136	do k=2,llm-1
1137	omdn=max(0.0,omega(k+1))
1138	omup=min(0.0,omega( k ))
1139
1140	! d_u_va(k) = -omdndudp(k)-omupdudp(k-1)
1141	! d_v_va(k) = -omdndvdp(k)-omupdvdp(k-1)
1142	d_q_va(k,1)= -omdndqdp(k)-omupdqdp(k-1)
1143	d_t_va(k) = -omdndtdp(k)-omupdtdp(k-1)+(omup+omdn)*cor(k)
1144	enddo
1145
1146	omdn=max(0.0,omega( 2 ))
1147	omup=min(0.0,omega(llm))
1148	! d_u_va( 1 ) = -omdn*dudp( 1 )
1149	! d_u_va(llm) = -omup*dudp(llm)
1150	! d_v_va( 1 ) = -omdn*dvdp( 1 )
1151	! d_v_va(llm) = -omup*dvdp(llm)
1152	d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1153	d_q_va(llm,1) = -omup*dqdp(llm)
1154	d_t_va( 1 ) = -omdndtdp( 1 )+omdncor( 1 )
1155	d_t_va(llm) = -omupdtdp(llm)!+omupcor(llm)
1156
1157	! if(abs(rlat(1))>10.) then
1158	! Calculate the tendency due agestrophic motions
1159	! du_age = fcoriolis*(v-vg)
1160	! dv_age = fcoriolis*(ug-u)
1161	! endif
1162
1163	! call writefield_phy('d_t_va',d_t_va,llm)
1164
1165	return
1166	end
1167
1168	!======================================================================
1169	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &
1170	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
1171	& ,ht_toga,vt_toga,hq_toga,vq_toga)
1172	implicit none
1173
1174	!-------------------------------------------------------------------------
1175	! Read TOGA-COARE forcing data
1176	!-------------------------------------------------------------------------
1177
1178	integer nlev_toga,nt_toga
1179	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1180	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1181	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1182	real w_toga(nlev_toga,nt_toga)
1183	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1184	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1185	character*80 fich_toga
1186
1187	integer k,ip
1188	real bid
1189
1190	integer iy,im,id,ih
1191
1192	real plev_min
1193
1194	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1195
1196	open(21,file=trim(fich_toga),form='formatted')
1197	read(21,'(a)')
1198	do ip = 1, nt_toga
1199	read(21,'(a)')
1200	read(21,'(a)')
1201	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1202	read(21,'(a)')
1203	read(21,'(a)')
1204
1205	do k = 1, nlev_toga
1206	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &
1207	& ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &
1208	& ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1209
1210	! conversion in SI units:
1211	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
1212	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
1213	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1214	! no water vapour tendency above 55 hPa
1215	if (plev_toga(k,ip) .lt. plev_min) then
1216	q_toga(k,ip) = 0.
1217	hq_toga(k,ip) = 0.
1218	vq_toga(k,ip) =0.
1219	endif
1220	enddo
1221
1222	ts_toga(ip)=ts_toga(ip)+273.15 ! K
1223	enddo
1224	close(21)
1225
1226	223 format(4i3,6f8.2)
1227	230 format(6f9.3,4e11.3)
1228
1229	return
1230	end
1231
1232	!-------------------------------------------------------------------------
1233	SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
1234	implicit none
1235
1236	!-------------------------------------------------------------------------
1237	! Read I.SANDU case forcing data
1238	!-------------------------------------------------------------------------
1239
1240	integer nlev_sandu,nt_sandu
1241	real ts_sandu(nt_sandu)
1242	character*80 fich_sandu
1243
1244	integer ip
1245	integer iy,im,id,ih
1246
1247	real plev_min
1248
1249	print*,'nlev_sandu',nlev_sandu
1250	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1251
1252	open(21,file=trim(fich_sandu),form='formatted')
1253	read(21,'(a)')
1254	do ip = 1, nt_sandu
1255	read(21,'(a)')
1256	read(21,'(a)')
1257	read(21,223) iy, im, id, ih, ts_sandu(ip)
1258	print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
1259	enddo
1260	close(21)
1261
1262	223 format(4i3,f8.2)
1263
1264	return
1265	end
1266
1267	!=====================================================================
1268	!-------------------------------------------------------------------------
1269	SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, &
1270	& ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
1271	implicit none
1272
1273	!-------------------------------------------------------------------------
1274	! Read Astex case forcing data
1275	!-------------------------------------------------------------------------
1276
1277	integer nlev_astex,nt_astex
1278	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
1279	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
1280	character*80 fich_astex
1281
1282	integer ip
1283	integer iy,im,id,ih
1284
1285	real plev_min
1286
1287	print*,'nlev_astex',nlev_astex
1288	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1289
1290	open(21,file=trim(fich_astex),form='formatted')
1291	read(21,'(a)')
1292	read(21,'(a)')
1293	do ip = 1, nt_astex
1294	read(21,'(a)')
1295	read(21,'(a)')
1296	read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), &
1297	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
1298	ts_astex(ip)=ts_astex(ip)+273.15
1299	print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), &
1300	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
1301	enddo
1302	close(21)
1303
1304	223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
1305
1306	return
1307	end
1308	!=====================================================================
1309	subroutine read_twpice(fich_twpice,nlevel,ntime &
1310	& ,T_srf,plev,T,q,u,v,omega &
1311	& ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1312
1313	!program reading forcings of the TWP-ICE experiment
1314
1315	! use netcdf
1316
1317	implicit none
1318
1319	#include "netcdf.inc"
1320
1321	integer ntime,nlevel
1322	integer l,k
1323	character*80 :: fich_twpice
1324	real*8 time(ntime)
1325	real*8 lat, lon, alt, phis
1326	real*8 lev(nlevel)
1327	real*8 plev(nlevel,ntime)
1328
1329	real*8 T(nlevel,ntime)
1330	real*8 q(nlevel,ntime),u(nlevel,ntime)
1331	real*8 v(nlevel,ntime)
1332	real*8 omega(nlevel,ntime), div(nlevel,ntime)
1333	real*8 T_adv_h(nlevel,ntime)
1334	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1335	real*8 q_adv_v(nlevel,ntime)
1336	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1337	real*8 s_adv_v(nlevel,ntime)
1338	real*8 p_srf_aver(ntime), p_srf_center(ntime)
1339	real*8 T_srf(ntime)
1340
1341	integer nid, ierr
1342	integer nbvar3d
1343	parameter(nbvar3d=20)
1344	integer var3didin(nbvar3d)
1345
1346	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1347	if (ierr.NE.NF_NOERR) then
1348	write(,) 'ERROR: Pb opening forcings cdf file '
1349	write(,) NF_STRERROR(ierr)
1350	stop ""
1351	endif
1352
1353	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1354	if(ierr/=NF_NOERR) then
1355	write(,) NF_STRERROR(ierr)
1356	stop 'lat'
1357	endif
1358
1359	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1360	if(ierr/=NF_NOERR) then
1361	write(,) NF_STRERROR(ierr)
1362	stop 'lon'
1363	endif
1364
1365	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1366	if(ierr/=NF_NOERR) then
1367	write(,) NF_STRERROR(ierr)
1368	stop 'alt'
1369	endif
1370
1371	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1372	if(ierr/=NF_NOERR) then
1373	write(,) NF_STRERROR(ierr)
1374	stop 'phis'
1375	endif
1376
1377	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1378	if(ierr/=NF_NOERR) then
1379	write(,) NF_STRERROR(ierr)
1380	stop 'T'
1381	endif
1382
1383	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1384	if(ierr/=NF_NOERR) then
1385	write(,) NF_STRERROR(ierr)
1386	stop 'q'
1387	endif
1388
1389	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1390	if(ierr/=NF_NOERR) then
1391	write(,) NF_STRERROR(ierr)
1392	stop 'u'
1393	endif
1394
1395	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1396	if(ierr/=NF_NOERR) then
1397	write(,) NF_STRERROR(ierr)
1398	stop 'v'
1399	endif
1400
1401	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1402	if(ierr/=NF_NOERR) then
1403	write(,) NF_STRERROR(ierr)
1404	stop 'omega'
1405	endif
1406
1407	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1408	if(ierr/=NF_NOERR) then
1409	write(,) NF_STRERROR(ierr)
1410	stop 'div'
1411	endif
1412
1413	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1414	if(ierr/=NF_NOERR) then
1415	write(,) NF_STRERROR(ierr)
1416	stop 'T_adv_h'
1417	endif
1418
1419	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1420	if(ierr/=NF_NOERR) then
1421	write(,) NF_STRERROR(ierr)
1422	stop 'T_adv_v'
1423	endif
1424
1425	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1426	if(ierr/=NF_NOERR) then
1427	write(,) NF_STRERROR(ierr)
1428	stop 'q_adv_h'
1429	endif
1430
1431	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1432	if(ierr/=NF_NOERR) then
1433	write(,) NF_STRERROR(ierr)
1434	stop 'q_adv_v'
1435	endif
1436
1437	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1438	if(ierr/=NF_NOERR) then
1439	write(,) NF_STRERROR(ierr)
1440	stop 's'
1441	endif
1442
1443	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1444	if(ierr/=NF_NOERR) then
1445	write(,) NF_STRERROR(ierr)
1446	stop 's_adv_h'
1447	endif
1448
1449	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1450	if(ierr/=NF_NOERR) then
1451	write(,) NF_STRERROR(ierr)
1452	stop 's_adv_v'
1453	endif
1454
1455	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1456	if(ierr/=NF_NOERR) then
1457	write(,) NF_STRERROR(ierr)
1458	stop 'p_srf_aver'
1459	endif
1460
1461	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1462	if(ierr/=NF_NOERR) then
1463	write(,) NF_STRERROR(ierr)
1464	stop 'p_srf_center'
1465	endif
1466
1467	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1468	if(ierr/=NF_NOERR) then
1469	write(,) NF_STRERROR(ierr)
1470	stop 'T_srf'
1471	endif
1472
1473	!dimensions lecture
1474	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1475
1476	!pressure
1477	do l=1,ntime
1478	do k=1,nlevel
1479	plev(k,l)=lev(k)
1480	enddo
1481	enddo
1482
1483	#ifdef NC_DOUBLE
1484	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1485	#else
1486	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1487	#endif
1488	if(ierr/=NF_NOERR) then
1489	write(,) NF_STRERROR(ierr)
1490	stop "getvarup"
1491	endif
1492	! write(,)'lecture lat ok',lat
1493
1494	#ifdef NC_DOUBLE
1495	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1496	#else
1497	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1498	#endif
1499	if(ierr/=NF_NOERR) then
1500	write(,) NF_STRERROR(ierr)
1501	stop "getvarup"
1502	endif
1503	! write(,)'lecture lon ok',lon
1504
1505	#ifdef NC_DOUBLE
1506	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1507	#else
1508	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1509	#endif
1510	if(ierr/=NF_NOERR) then
1511	write(,) NF_STRERROR(ierr)
1512	stop "getvarup"
1513	endif
1514	! write(,)'lecture alt ok',alt
1515
1516	#ifdef NC_DOUBLE
1517	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1518	#else
1519	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1520	#endif
1521	if(ierr/=NF_NOERR) then
1522	write(,) NF_STRERROR(ierr)
1523	stop "getvarup"
1524	endif
1525	! write(,)'lecture phis ok',phis
1526
1527	#ifdef NC_DOUBLE
1528	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1529	#else
1530	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1531	#endif
1532	if(ierr/=NF_NOERR) then
1533	write(,) NF_STRERROR(ierr)
1534	stop "getvarup"
1535	endif
1536	! write(,)'lecture T ok'
1537
1538	#ifdef NC_DOUBLE
1539	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1540	#else
1541	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1542	#endif
1543	if(ierr/=NF_NOERR) then
1544	write(,) NF_STRERROR(ierr)
1545	stop "getvarup"
1546	endif
1547	! write(,)'lecture q ok'
1548	!q in kg/kg
1549	do l=1,ntime
1550	do k=1,nlevel
1551	q(k,l)=q(k,l)/1000.
1552	enddo
1553	enddo
1554	#ifdef NC_DOUBLE
1555	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1556	#else
1557	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1558	#endif
1559	if(ierr/=NF_NOERR) then
1560	write(,) NF_STRERROR(ierr)
1561	stop "getvarup"
1562	endif
1563	! write(,)'lecture u ok'
1564
1565	#ifdef NC_DOUBLE
1566	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1567	#else
1568	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1569	#endif
1570	if(ierr/=NF_NOERR) then
1571	write(,) NF_STRERROR(ierr)
1572	stop "getvarup"
1573	endif
1574	! write(,)'lecture v ok'
1575
1576	#ifdef NC_DOUBLE
1577	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1578	#else
1579	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1580	#endif
1581	if(ierr/=NF_NOERR) then
1582	write(,) NF_STRERROR(ierr)
1583	stop "getvarup"
1584	endif
1585	! write(,)'lecture omega ok'
1586	!omega in mb/hour
1587	do l=1,ntime
1588	do k=1,nlevel
1589	omega(k,l)=omega(k,l)*100./3600.
1590	enddo
1591	enddo
1592
1593	#ifdef NC_DOUBLE
1594	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1595	#else
1596	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1597	#endif
1598	if(ierr/=NF_NOERR) then
1599	write(,) NF_STRERROR(ierr)
1600	stop "getvarup"
1601	endif
1602	! write(,)'lecture div ok'
1603
1604	#ifdef NC_DOUBLE
1605	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1606	#else
1607	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1608	#endif
1609	if(ierr/=NF_NOERR) then
1610	write(,) NF_STRERROR(ierr)
1611	stop "getvarup"
1612	endif
1613	! write(,)'lecture T_adv_h ok'
1614	!T adv in K/s
1615	do l=1,ntime
1616	do k=1,nlevel
1617	T_adv_h(k,l)=T_adv_h(k,l)/3600.
1618	enddo
1619	enddo
1620
1621
1622	#ifdef NC_DOUBLE
1623	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1624	#else
1625	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1626	#endif
1627	if(ierr/=NF_NOERR) then
1628	write(,) NF_STRERROR(ierr)
1629	stop "getvarup"
1630	endif
1631	! write(,)'lecture T_adv_v ok'
1632	!T adv in K/s
1633	do l=1,ntime
1634	do k=1,nlevel
1635	T_adv_v(k,l)=T_adv_v(k,l)/3600.
1636	enddo
1637	enddo
1638
1639	#ifdef NC_DOUBLE
1640	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1641	#else
1642	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1643	#endif
1644	if(ierr/=NF_NOERR) then
1645	write(,) NF_STRERROR(ierr)
1646	stop "getvarup"
1647	endif
1648	! write(,)'lecture q_adv_h ok'
1649	!q adv in kg/kg/s
1650	do l=1,ntime
1651	do k=1,nlevel
1652	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1653	enddo
1654	enddo
1655
1656
1657	#ifdef NC_DOUBLE
1658	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1659	#else
1660	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1661	#endif
1662	if(ierr/=NF_NOERR) then
1663	write(,) NF_STRERROR(ierr)
1664	stop "getvarup"
1665	endif
1666	! write(,)'lecture q_adv_v ok'
1667	!q adv in kg/kg/s
1668	do l=1,ntime
1669	do k=1,nlevel
1670	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1671	enddo
1672	enddo
1673
1674
1675	#ifdef NC_DOUBLE
1676	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1677	#else
1678	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1679	#endif
1680	if(ierr/=NF_NOERR) then
1681	write(,) NF_STRERROR(ierr)
1682	stop "getvarup"
1683	endif
1684
1685	#ifdef NC_DOUBLE
1686	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1687	#else
1688	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1689	#endif
1690	if(ierr/=NF_NOERR) then
1691	write(,) NF_STRERROR(ierr)
1692	stop "getvarup"
1693	endif
1694
1695	#ifdef NC_DOUBLE
1696	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1697	#else
1698	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1699	#endif
1700	if(ierr/=NF_NOERR) then
1701	write(,) NF_STRERROR(ierr)
1702	stop "getvarup"
1703	endif
1704
1705	#ifdef NC_DOUBLE
1706	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1707	#else
1708	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1709	#endif
1710	if(ierr/=NF_NOERR) then
1711	write(,) NF_STRERROR(ierr)
1712	stop "getvarup"
1713	endif
1714
1715	#ifdef NC_DOUBLE
1716	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1717	#else
1718	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1719	#endif
1720	if(ierr/=NF_NOERR) then
1721	write(,) NF_STRERROR(ierr)
1722	stop "getvarup"
1723	endif
1724
1725	#ifdef NC_DOUBLE
1726	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1727	#else
1728	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
1729	#endif
1730	if(ierr/=NF_NOERR) then
1731	write(,) NF_STRERROR(ierr)
1732	stop "getvarup"
1733	endif
1734	! write(,)'lecture T_srf ok', T_srf
1735
1736	return
1737	end subroutine read_twpice
1738	!=====================================================================
1739	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
1740
1741	! use netcdf
1742
1743	implicit none
1744	#include "netcdf.inc"
1745	integer nid,ttm,llm
1746	real*8 time(ttm)
1747	real*8 lev(llm)
1748	integer ierr
1749
1750	integer timevar,levvar
1751	integer timelen,levlen
1752	integer timedimin,levdimin
1753
1754	! Control & lecture on dimensions
1755	! ===============================
1756	ierr=NF_INQ_DIMID(nid,"time",timedimin)
1757	ierr=NF_INQ_VARID(nid,"time",timevar)
1758	if (ierr.NE.NF_NOERR) then
1759	write(,) 'ERROR: Field <time> is missing'
1760	stop ""
1761	endif
1762	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
1763
1764	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
1765	ierr=NF_INQ_VARID(nid,"lev",levvar)
1766	if (ierr.NE.NF_NOERR) then
1767	write(,) 'ERROR: Field <lev> is lacking'
1768	stop ""
1769	endif
1770	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
1771
1772	if((timelen/=ttm).or.(levlen/=llm)) then
1773	write(,) 'ERROR: Not the good lenght for axis'
1774	write(,) 'longitude: ',timelen,ttm+1
1775	write(,) 'latitude: ',levlen,llm
1776	stop ""
1777	endif
1778
1779	!#ifdef NC_DOUBLE
1780	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
1781	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
1782	!#else
1783	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
1784	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
1785	!#endif
1786
1787	return
1788	end
1789	!=====================================================================
1790
1791	SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof &
1792	& ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof &
1793	& ,omega_prof,o3mmr_prof &
1794	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
1795	& ,omega_mod,o3mmr_mod,mxcalc)
1796
1797	implicit none
1798
1799	#include "dimensions.h"
1800
1801	!-------------------------------------------------------------------------
1802	! Vertical interpolation of SANDUREF forcing data onto model levels
1803	!-------------------------------------------------------------------------
1804
1805	integer nlevmax
1806	parameter (nlevmax=41)
1807	integer nlev_sandu,mxcalc
1808	! real play(llm), plev_prof(nlevmax)
1809	! real t_prof(nlevmax),q_prof(nlevmax)
1810	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1811	! real ht_prof(nlevmax),vt_prof(nlevmax)
1812	! real hq_prof(nlevmax),vq_prof(nlevmax)
1813
1814	real play(llm), plev_prof(nlev_sandu)
1815	real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
1816	real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
1817	real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
1818
1819	real t_mod(llm),thl_mod(llm),q_mod(llm)
1820	real u_mod(llm),v_mod(llm), w_mod(llm)
1821	real omega_mod(llm),o3mmr_mod(llm)
1822
1823	integer l,k,k1,k2
1824	real frac,frac1,frac2,fact
1825
1826	do l = 1, llm
1827
1828	if (play(l).ge.plev_prof(nlev_sandu)) then
1829
1830	mxcalc=l
1831	k1=0
1832	k2=0
1833
1834	if (play(l).le.plev_prof(1)) then
1835
1836	do k = 1, nlev_sandu-1
1837	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
1838	k1=k
1839	k2=k+1
1840	endif
1841	enddo
1842
1843	if (k1.eq.0 .or. k2.eq.0) then
1844	write(,) 'PB! k1, k2 = ',k1,k2
1845	write(,) 'l,play(l) = ',l,play(l)/100
1846	do k = 1, nlev_sandu-1
1847	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1848	enddo
1849	endif
1850
1851	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1852	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
1853	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
1854	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
1855	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1856	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1857	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1858	omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
1859	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
1860
1861	else !play>plev_prof(1)
1862
1863	k1=1
1864	k2=2
1865	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1866	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1867	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
1868	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
1869	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
1870	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
1871	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
1872	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
1873	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
1874	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
1875
1876	endif ! play.le.plev_prof(1)
1877
1878	else ! above max altitude of forcing file
1879
1880	!jyg
1881	fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
1882	fact = max(fact,0.) !jyg
1883	fact = exp(-fact) !jyg
1884	t_mod(l)= t_prof(nlev_sandu) !jyg
1885	thl_mod(l)= thl_prof(nlev_sandu) !jyg
1886	q_mod(l)= q_prof(nlev_sandu)*fact !jyg
1887	u_mod(l)= u_prof(nlev_sandu)*fact !jyg
1888	v_mod(l)= v_prof(nlev_sandu)*fact !jyg
1889	w_mod(l)= w_prof(nlev_sandu)*fact !jyg
1890	omega_mod(l)= omega_prof(nlev_sandu)*fact !jyg
1891	o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact !jyg
1892
1893	endif ! play
1894
1895	enddo ! l
1896
1897	do l = 1,llm
1898	! print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
1899	! $ l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
1900	enddo
1901
1902	return
1903	end
1904	!=====================================================================
1905	SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof &
1906	& ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof &
1907	& ,w_prof,tke_prof,o3mmr_prof &
1908	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
1909	& ,tke_mod,o3mmr_mod,mxcalc)
1910
1911	implicit none
1912
1913	#include "dimensions.h"
1914
1915	!-------------------------------------------------------------------------
1916	! Vertical interpolation of Astex forcing data onto model levels
1917	!-------------------------------------------------------------------------
1918
1919	integer nlevmax
1920	parameter (nlevmax=41)
1921	integer nlev_astex,mxcalc
1922	! real play(llm), plev_prof(nlevmax)
1923	! real t_prof(nlevmax),qv_prof(nlevmax)
1924	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1925	! real ht_prof(nlevmax),vt_prof(nlevmax)
1926	! real hq_prof(nlevmax),vq_prof(nlevmax)
1927
1928	real play(llm), plev_prof(nlev_astex)
1929	real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
1930	real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
1931	real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
1932	real qt_prof(nlev_astex),tke_prof(nlev_astex)
1933
1934	real t_mod(llm),thl_mod(llm),qv_mod(llm)
1935	real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
1936	real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
1937
1938	integer l,k,k1,k2
1939	real frac,frac1,frac2,fact
1940
1941	do l = 1, llm
1942
1943	if (play(l).ge.plev_prof(nlev_astex)) then
1944
1945	mxcalc=l
1946	k1=0
1947	k2=0
1948
1949	if (play(l).le.plev_prof(1)) then
1950
1951	do k = 1, nlev_astex-1
1952	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
1953	k1=k
1954	k2=k+1
1955	endif
1956	enddo
1957
1958	if (k1.eq.0 .or. k2.eq.0) then
1959	write(,) 'PB! k1, k2 = ',k1,k2
1960	write(,) 'l,play(l) = ',l,play(l)/100
1961	do k = 1, nlev_astex-1
1962	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1963	enddo
1964	endif
1965
1966	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1967	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
1968	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
1969	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
1970	ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
1971	qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
1972	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1973	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1974	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1975	tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
1976	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
1977
1978	else !play>plev_prof(1)
1979
1980	k1=1
1981	k2=2
1982	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1983	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1984	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
1985	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
1986	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
1987	ql_mod(l)= frac1ql_prof(k1) - frac2ql_prof(k2)
1988	qt_mod(l)= frac1qt_prof(k1) - frac2qt_prof(k2)
1989	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
1990	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
1991	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
1992	tke_mod(l)= frac1tke_prof(k1) - frac2tke_prof(k2)
1993	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
1994
1995	endif ! play.le.plev_prof(1)
1996
1997	else ! above max altitude of forcing file
1998
1999	!jyg
2000	fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
2001	fact = max(fact,0.) !jyg
2002	fact = exp(-fact) !jyg
2003	t_mod(l)= t_prof(nlev_astex) !jyg
2004	thl_mod(l)= thl_prof(nlev_astex) !jyg
2005	qv_mod(l)= qv_prof(nlev_astex)*fact !jyg
2006	ql_mod(l)= ql_prof(nlev_astex)*fact !jyg
2007	qt_mod(l)= qt_prof(nlev_astex)*fact !jyg
2008	u_mod(l)= u_prof(nlev_astex)*fact !jyg
2009	v_mod(l)= v_prof(nlev_astex)*fact !jyg
2010	w_mod(l)= w_prof(nlev_astex)*fact !jyg
2011	tke_mod(l)= tke_prof(nlev_astex)*fact !jyg
2012	o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact !jyg
2013
2014	endif ! play
2015
2016	enddo ! l
2017
2018	do l = 1,llm
2019	! print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
2020	! $ l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
2021	enddo
2022
2023	return
2024	end
2025
2026	!======================================================================
2027	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play &
2028	& ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &
2029	& ,dth_dyn,dqh_dyn)
2030	implicit none
2031
2032	!-------------------------------------------------------------------------
2033	! Read RICO forcing data
2034	!-------------------------------------------------------------------------
2035	#include "dimensions.h"
2036
2037
2038	integer nlev_rico
2039	real ts_rico,ps_rico
2040	real t_rico(llm),q_rico(llm)
2041	real u_rico(llm),v_rico(llm)
2042	real w_rico(llm)
2043	real dth_dyn(llm)
2044	real dqh_dyn(llm)
2045
2046
2047	real play(llm),zlay(llm)
2048
2049
2050	real prico(nlev_rico),zrico(nlev_rico)
2051
2052	character*80 fich_rico
2053
2054	integer k,l
2055
2056
2057	print*,fich_rico
2058	open(21,file=trim(fich_rico),form='formatted')
2059	do k=1,llm
2060	zlay(k)=0.
2061	enddo
2062
2063	read(21,*) ps_rico,ts_rico
2064	prico(1)=ps_rico
2065	zrico(1)=0.0
2066	do l=2,nlev_rico
2067	read(21,*) k,prico(l),zrico(l)
2068	enddo
2069	close(21)
2070
2071	do k=1,llm
2072	do l=1,80
2073	if(prico(l)>play(k)) then
2074	if(play(k)>prico(l+1)) then
2075	zlay(k)=zrico(l)+(play(k)-prico(l)) * &
2076	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2077	else
2078	zlay(k)=zrico(l)+(play(k)-prico(80))* &
2079	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
2080	endif
2081	endif
2082	enddo
2083	print*,k,zlay(k)
2084	! U
2085	if(0 < zlay(k) .and. zlay(k) < 4000) then
2086	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2087	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2088	u_rico(k)= -1.9 + (30.0 + 1.9) / &
2089	& (12000 - 4000) * (zlay(k) - 4000)
2090	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2091	u_rico(k)=30.0
2092	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2093	u_rico(k)=30.0 - (30.0) / &
2094	& (20000 - 13000) * (zlay(k) - 13000)
2095	else
2096	u_rico(k)=0.0
2097	endif
2098
2099	!Q_v
2100	if(0 < zlay(k) .and. zlay(k) < 740) then
2101	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2102	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2103	q_rico(k)=13.8 + (2.4 - 13.8) / &
2104	& (3260 - 740) * (zlay(k) - 740)
2105	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2106	q_rico(k)=2.4 + (1.8 - 2.4) / &
2107	& (4000 - 3260) * (zlay(k) - 3260)
2108	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2109	q_rico(k)=1.8 + (0 - 1.8) / &
2110	& (10000 - 4000) * (zlay(k) - 4000)
2111	else
2112	q_rico(k)=0.0
2113	endif
2114
2115	!T
2116	if(0 < zlay(k) .and. zlay(k) < 740) then
2117	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2118	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2119	t_rico(k)=292.0 + (278.0 - 292.0) / &
2120	& (4000 - 740) * (zlay(k) - 740)
2121	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2122	t_rico(k)=278.0 + (203.0 - 278.0) / &
2123	& (15000 - 4000) * (zlay(k) - 4000)
2124	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2125	t_rico(k)=203.0 + (194.0 - 203.0) / &
2126	& (17500 - 15000)* (zlay(k) - 15000)
2127	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2128	t_rico(k)=194.0 + (206.0 - 194.0) / &
2129	& (20000 - 17500)* (zlay(k) - 17500)
2130	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2131	t_rico(k)=206.0 + (270.0 - 206.0) / &
2132	& (60000 - 20000)* (zlay(k) - 20000)
2133	endif
2134
2135	! W
2136	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2137	w_rico(k)=- (0.005/2260) * zlay(k)
2138	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2139	w_rico(k)=- 0.005
2140	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2141	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2142	else
2143	w_rico(k)=0.0
2144	endif
2145
2146	! dThrz+dTsw0+dTlw0
2147	if(0 < zlay(k) .and. zlay(k) < 4000) then
2148	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ &
2149	& (864004000) zlay(k)
2150	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2151	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) / &
2152	& (86400(5000 - 4000)) (zlay(k) - 4000)
2153	else
2154	dth_dyn(k)=0.0
2155	endif
2156	! dQhrz
2157	if(0 < zlay(k) .and. zlay(k) < 3000) then
2158	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/ &
2159	& (864003000) (zlay(k))
2160	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2161	dqh_dyn(k)=0.345 / 86400
2162	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2163	dqh_dyn(k)=0.345 / 86400 + &
2164	& (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2165	else
2166	dqh_dyn(k)=0.0
2167	endif
2168
2169	!? if(play(k)>6e4) then
2170	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2171	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
2172	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2173	!? *(6e4-play(k))/(6e4-3e4)
2174	!? else
2175	!? ratqs0(1,k)=ratqshaut
2176	!? endif
2177
2178	enddo
2179
2180	do k=1,llm
2181	q_rico(k)=q_rico(k)/1e3
2182	dqh_dyn(k)=dqh_dyn(k)/1e3
2183	v_rico(k)=-3.8
2184	enddo
2185
2186	return
2187	end
2188
2189	!======================================================================
2190	SUBROUTINE interp_sandu_time(day,day1,annee_ref &
2191	& ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu &
2192	& ,nlev_sandu,ts_sandu,ts_prof)
2193	implicit none
2194
2195	!---------------------------------------------------------------------------------------
2196	! Time interpolation of a 2D field to the timestep corresponding to day
2197	!
2198	! day: current julian day (e.g. 717538.2)
2199	! day1: first day of the simulation
2200	! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
2201	! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
2202	!---------------------------------------------------------------------------------------
2203	! inputs:
2204	integer annee_ref
2205	integer nt_sandu,nlev_sandu
2206	integer year_ini_sandu
2207	real day, day1,day_ini_sandu,dt_sandu
2208	real ts_sandu(nt_sandu)
2209	! outputs:
2210	real ts_prof
2211	! local:
2212	integer it_sandu1, it_sandu2
2213	real timeit,time_sandu1,time_sandu2,frac
2214	! Check that initial day of the simulation consistent with SANDU period:
2215	if (annee_ref.ne.2006 ) then
2216	print*,'Pour SANDUREF, annee_ref doit etre 2006 '
2217	print*,'Changer annee_ref dans run.def'
2218	stop
2219	endif
2220	! if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
2221	! print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
2222	! print*,'Changer dayref dans run.def'
2223	! stop
2224	! endif
2225
2226	! Determine timestep relative to the 1st day of TOGA-COARE:
2227	! timeit=(day-day1)*86400.
2228	! if (annee_ref.eq.1992) then
2229	! timeit=(day-day_ini_sandu)*86400.
2230	! else
2231	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2232	! endif
2233	timeit=(day-day_ini_sandu)*86400
2234
2235	! Determine the closest observation times:
2236	it_sandu1=INT(timeit/dt_sandu)+1
2237	it_sandu2=it_sandu1 + 1
2238	time_sandu1=(it_sandu1-1)*dt_sandu
2239	time_sandu2=(it_sandu2-1)*dt_sandu
2240	print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
2241	print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2', &
2242	& it_sandu1,it_sandu2,time_sandu1,time_sandu2
2243
2244	if (it_sandu1 .ge. nt_sandu) then
2245	write(,) 'PB-stop: day, it_sandu1, it_sandu2, timeit: ' &
2246	& ,day,it_sandu1,it_sandu2,timeit/86400.
2247	stop
2248	endif
2249
2250	! time interpolation:
2251	frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
2252	frac=max(frac,0.0)
2253
2254	ts_prof = ts_sandu(it_sandu2) &
2255	& -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
2256
2257	print*, &
2258	&'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:', &
2259	&day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1, &
2260	&it_sandu2,ts_prof
2261
2262	return
2263	END
2264	!=====================================================================
2265	!-------------------------------------------------------------------------
2266	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu, &
2267	& sens,flat,adv_theta,rad_theta,adv_qt)
2268	implicit none
2269
2270	!-------------------------------------------------------------------------
2271	! Read ARM_CU case forcing data
2272	!-------------------------------------------------------------------------
2273
2274	integer nlev_armcu,nt_armcu
2275	real sens(nt_armcu),flat(nt_armcu)
2276	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2277	character*80 fich_armcu
2278
2279	integer ip
2280
2281	integer iy,im,id,ih,in
2282
2283	print*,'nlev_armcu',nlev_armcu
2284
2285	open(21,file=trim(fich_armcu),form='formatted')
2286	read(21,'(a)')
2287	do ip = 1, nt_armcu
2288	read(21,'(a)')
2289	read(21,'(a)')
2290	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip), &
2291	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
2292	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip), &
2293	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
2294	enddo
2295	close(21)
2296
2297	223 format(5i3,5f8.3)
2298
2299	return
2300	end
2301
2302	!=====================================================================
2303	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof &
2304	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
2305	& ,ht_prof,vt_prof,hq_prof,vq_prof &
2306	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
2307	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2308
2309	implicit none
2310
2311	#include "dimensions.h"
2312
2313	!-------------------------------------------------------------------------
2314	! Vertical interpolation of TOGA-COARE forcing data onto model levels
2315	!-------------------------------------------------------------------------
2316
2317	integer nlevmax
2318	parameter (nlevmax=41)
2319	integer nlev_toga,mxcalc
2320	! real play(llm), plev_prof(nlevmax)
2321	! real t_prof(nlevmax),q_prof(nlevmax)
2322	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2323	! real ht_prof(nlevmax),vt_prof(nlevmax)
2324	! real hq_prof(nlevmax),vq_prof(nlevmax)
2325
2326	real play(llm), plev_prof(nlev_toga)
2327	real t_prof(nlev_toga),q_prof(nlev_toga)
2328	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2329	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2330	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2331
2332	real t_mod(llm),q_mod(llm)
2333	real u_mod(llm),v_mod(llm), w_mod(llm)
2334	real ht_mod(llm),vt_mod(llm)
2335	real hq_mod(llm),vq_mod(llm)
2336
2337	integer l,k,k1,k2
2338	real frac,frac1,frac2,fact
2339
2340	do l = 1, llm
2341
2342	if (play(l).ge.plev_prof(nlev_toga)) then
2343
2344	mxcalc=l
2345	k1=0
2346	k2=0
2347
2348	if (play(l).le.plev_prof(1)) then
2349
2350	do k = 1, nlev_toga-1
2351	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2352	k1=k
2353	k2=k+1
2354	endif
2355	enddo
2356
2357	if (k1.eq.0 .or. k2.eq.0) then
2358	write(,) 'PB! k1, k2 = ',k1,k2
2359	write(,) 'l,play(l) = ',l,play(l)/100
2360	do k = 1, nlev_toga-1
2361	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2362	enddo
2363	endif
2364
2365	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2366	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2367	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2368	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2369	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2370	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2371	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2372	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2373	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2374	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2375
2376	else !play>plev_prof(1)
2377
2378	k1=1
2379	k2=2
2380	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2381	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2382	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2383	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
2384	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2385	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2386	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2387	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
2388	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
2389	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
2390	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
2391
2392	endif ! play.le.plev_prof(1)
2393
2394	else ! above max altitude of forcing file
2395
2396	!jyg
2397	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2398	fact = max(fact,0.) !jyg
2399	fact = exp(-fact) !jyg
2400	t_mod(l)= t_prof(nlev_toga) !jyg
2401	q_mod(l)= q_prof(nlev_toga)*fact !jyg
2402	u_mod(l)= u_prof(nlev_toga)*fact !jyg
2403	v_mod(l)= v_prof(nlev_toga)*fact !jyg
2404	w_mod(l)= 0.0 !jyg
2405	ht_mod(l)= ht_prof(nlev_toga) !jyg
2406	vt_mod(l)= vt_prof(nlev_toga) !jyg
2407	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
2408	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
2409
2410	endif ! play
2411
2412	enddo ! l
2413
2414	! do l = 1,llm
2415	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2416	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2417	! enddo
2418
2419	return
2420	end
2421
2422	!======================================================================
2423	SUBROUTINE interp_astex_time(day,day1,annee_ref &
2424	& ,year_ini_astex,day_ini_astex,nt_astex,dt_astex &
2425	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
2426	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
2427	& ,ufa_prof,vfa_prof)
2428	implicit none
2429
2430	!---------------------------------------------------------------------------------------
2431	! Time interpolation of a 2D field to the timestep corresponding to day
2432	!
2433	! day: current julian day (e.g. 717538.2)
2434	! day1: first day of the simulation
2435	! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
2436	! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
2437	!---------------------------------------------------------------------------------------
2438
2439	! inputs:
2440	integer annee_ref
2441	integer nt_astex,nlev_astex
2442	integer year_ini_astex
2443	real day, day1,day_ini_astex,dt_astex
2444	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
2445	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
2446	! outputs:
2447	real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
2448	! local:
2449	integer it_astex1, it_astex2
2450	real timeit,time_astex1,time_astex2,frac
2451
2452	! Check that initial day of the simulation consistent with ASTEX period:
2453	if (annee_ref.ne.1992 ) then
2454	print*,'Pour Astex, annee_ref doit etre 1992 '
2455	print*,'Changer annee_ref dans run.def'
2456	stop
2457	endif
2458	if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
2459	print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
2460	print*,'Changer dayref dans run.def'
2461	stop
2462	endif
2463
2464	! Determine timestep relative to the 1st day of TOGA-COARE:
2465	! timeit=(day-day1)*86400.
2466	! if (annee_ref.eq.1992) then
2467	! timeit=(day-day_ini_astex)*86400.
2468	! else
2469	! timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
2470	! endif
2471	timeit=(day-day_ini_astex)*86400
2472
2473	! Determine the closest observation times:
2474	it_astex1=INT(timeit/dt_astex)+1
2475	it_astex2=it_astex1 + 1
2476	time_astex1=(it_astex1-1)*dt_astex
2477	time_astex2=(it_astex2-1)*dt_astex
2478	print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
2479	print *,'it_astex1,it_astex2,time_astex1,time_astex2', &
2480	& it_astex1,it_astex2,time_astex1,time_astex2
2481
2482	if (it_astex1 .ge. nt_astex) then
2483	write(,) 'PB-stop: day, it_astex1, it_astex2, timeit: ' &
2484	& ,day,it_astex1,it_astex2,timeit/86400.
2485	stop
2486	endif
2487
2488	! time interpolation:
2489	frac=(time_astex2-timeit)/(time_astex2-time_astex1)
2490	frac=max(frac,0.0)
2491
2492	div_prof = div_astex(it_astex2) &
2493	& -frac*(div_astex(it_astex2)-div_astex(it_astex1))
2494	ts_prof = ts_astex(it_astex2) &
2495	& -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
2496	ug_prof = ug_astex(it_astex2) &
2497	& -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
2498	vg_prof = vg_astex(it_astex2) &
2499	& -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
2500	ufa_prof = ufa_astex(it_astex2) &
2501	& -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
2502	vfa_prof = vfa_astex(it_astex2) &
2503	& -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
2504
2505	print*, &
2506	&'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:', &
2507	&day,annee_ref,day_ini_astex,timeit/86400.,it_astex1, &
2508	&it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
2509
2510	return
2511	END
2512
2513	!======================================================================
2514	SUBROUTINE interp_toga_time(day,day1,annee_ref &
2515	& ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga &
2516	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
2517	& ,ht_toga,vt_toga,hq_toga,vq_toga &
2518	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
2519	& ,ht_prof,vt_prof,hq_prof,vq_prof)
2520	implicit none
2521
2522	!---------------------------------------------------------------------------------------
2523	! Time interpolation of a 2D field to the timestep corresponding to day
2524	!
2525	! day: current julian day (e.g. 717538.2)
2526	! day1: first day of the simulation
2527	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
2528	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
2529	!---------------------------------------------------------------------------------------
2530
2531	#include "compar1d.h"
2532
2533	! inputs:
2534	integer annee_ref
2535	integer nt_toga,nlev_toga
2536	integer year_ini_toga
2537	real day, day1,day_ini_toga,dt_toga
2538	real ts_toga(nt_toga)
2539	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
2540	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
2541	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
2542	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
2543	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
2544	! outputs:
2545	real ts_prof
2546	real plev_prof(nlev_toga),t_prof(nlev_toga)
2547	real q_prof(nlev_toga),u_prof(nlev_toga)
2548	real v_prof(nlev_toga),w_prof(nlev_toga)
2549	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2550	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2551	! local:
2552	integer it_toga1, it_toga2,k
2553	real timeit,time_toga1,time_toga2,frac
2554
2555
2556	if (forcing_type.eq.2) then
2557	! Check that initial day of the simulation consistent with TOGA-COARE period:
2558	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
2559	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
2560	print*,'Changer annee_ref dans run.def'
2561	stop
2562	endif
2563	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
2564	print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
2565	print*,'Changer dayref dans run.def'
2566	stop
2567	endif
2568	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
2569	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
2570	print*,'Changer dayref ou nday dans run.def'
2571	stop
2572	endif
2573
2574	else if (forcing_type.eq.4) then
2575
2576	! Check that initial day of the simulation consistent with TWP-ICE period:
2577	if (annee_ref.ne.2006) then
2578	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
2579	print*,'Changer annee_ref dans run.def'
2580	stop
2581	endif
2582	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
2583	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
2584	print*,'Changer dayref dans run.def'
2585	stop
2586	endif
2587	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
2588	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
2589	print*,'Changer dayref ou nday dans run.def'
2590	stop
2591	endif
2592
2593	endif
2594
2595	! Determine timestep relative to the 1st day of TOGA-COARE:
2596	! timeit=(day-day1)*86400.
2597	! if (annee_ref.eq.1992) then
2598	! timeit=(day-day_ini_toga)*86400.
2599	! else
2600	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2601	! endif
2602	timeit=(day-day_ini_toga)*86400
2603
2604	! Determine the closest observation times:
2605	it_toga1=INT(timeit/dt_toga)+1
2606	it_toga2=it_toga1 + 1
2607	time_toga1=(it_toga1-1)*dt_toga
2608	time_toga2=(it_toga2-1)*dt_toga
2609
2610	if (it_toga1 .ge. nt_toga) then
2611	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: ' &
2612	& ,day,it_toga1,it_toga2,timeit/86400.
2613	stop
2614	endif
2615
2616	! time interpolation:
2617	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
2618	frac=max(frac,0.0)
2619
2620	ts_prof = ts_toga(it_toga2) &
2621	& -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
2622
2623	! print*,
2624	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
2625	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
2626
2627	do k=1,nlev_toga
2628	plev_prof(k) = 100.*(plev_toga(k,it_toga2) &
2629	& -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
2630	t_prof(k) = t_toga(k,it_toga2) &
2631	& -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
2632	q_prof(k) = q_toga(k,it_toga2) &
2633	& -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
2634	u_prof(k) = u_toga(k,it_toga2) &
2635	& -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
2636	v_prof(k) = v_toga(k,it_toga2) &
2637	& -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
2638	w_prof(k) = w_toga(k,it_toga2) &
2639	& -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
2640	ht_prof(k) = ht_toga(k,it_toga2) &
2641	& -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
2642	vt_prof(k) = vt_toga(k,it_toga2) &
2643	& -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
2644	hq_prof(k) = hq_toga(k,it_toga2) &
2645	& -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
2646	vq_prof(k) = vq_toga(k,it_toga2) &
2647	& -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
2648	enddo
2649
2650	return
2651	END
2652
2653	!======================================================================
2654	SUBROUTINE interp_armcu_time(day,day1,annee_ref &
2655	& ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu &
2656	& ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu &
2657	& ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
2658	implicit none
2659
2660	!---------------------------------------------------------------------------------------
2661	! Time interpolation of a 2D field to the timestep corresponding to day
2662	!
2663	! day: current julian day (e.g. 717538.2)
2664	! day1: first day of the simulation
2665	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
2666	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
2667	! fs= sensible flux
2668	! fl= latent flux
2669	! at,rt,aqt= advective and radiative tendencies
2670	!---------------------------------------------------------------------------------------
2671
2672	! inputs:
2673	integer annee_ref
2674	integer nt_armcu,nlev_armcu
2675	integer year_ini_armcu
2676	real day, day1,day_ini_armcu,dt_armcu
2677	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
2678	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
2679	! outputs:
2680	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2681	! local:
2682	integer it_armcu1, it_armcu2,k
2683	real timeit,time_armcu1,time_armcu2,frac
2684
2685	! Check that initial day of the simulation consistent with ARMCU period:
2686	if (annee_ref.ne.1997 ) then
2687	print*,'Pour ARMCU, annee_ref doit etre 1997 '
2688	print*,'Changer annee_ref dans run.def'
2689	stop
2690	endif
2691
2692	timeit=(day-day_ini_armcu)*86400
2693
2694	! Determine the closest observation times:
2695	it_armcu1=INT(timeit/dt_armcu)+1
2696	it_armcu2=it_armcu1 + 1
2697	time_armcu1=(it_armcu1-1)*dt_armcu
2698	time_armcu2=(it_armcu2-1)*dt_armcu
2699	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
2700	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', &
2701	& it_armcu1,it_armcu2,time_armcu1,time_armcu2
2702
2703	if (it_armcu1 .ge. nt_armcu) then
2704	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' &
2705	& ,day,it_armcu1,it_armcu2,timeit/86400.
2706	stop
2707	endif
2708
2709	! time interpolation:
2710	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
2711	frac=max(frac,0.0)
2712
2713	fs_prof = fs_armcu(it_armcu2) &
2714	& -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
2715	fl_prof = fl_armcu(it_armcu2) &
2716	& -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
2717	at_prof = at_armcu(it_armcu2) &
2718	& -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
2719	rt_prof = rt_armcu(it_armcu2) &
2720	& -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
2721	aqt_prof = aqt_armcu(it_armcu2) &
2722	& -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
2723
2724	print*, &
2725	&'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', &
2726	&day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, &
2727	&it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2728
2729	return
2730	END
2731
2732	!=====================================================================
2733	subroutine readprofiles(nlev_max,kmax,ntrac,height, &
2734	& thlprof,qtprof,uprof, &
2735	& vprof,e12prof,ugprof,vgprof, &
2736	& wfls,dqtdxls,dqtdyls,dqtdtls, &
2737	& thlpcar,tracer,nt1,nt2)
2738	implicit none
2739
2740	integer nlev_max,kmax,kmax2,ntrac
2741	logical :: llesread = .true.
2742
2743	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), &
2744	& uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), &
2745	& ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), &
2746	& dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), &
2747	& thlpcar(nlev_max),tracer(nlev_max,ntrac)
2748
2749	integer, parameter :: ilesfile=1
2750	integer :: ierr,k,itrac,nt1,nt2
2751
2752	if(.not.(llesread)) return
2753
2754	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2755	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2756	read (ilesfile,*) kmax
2757	do k=1,kmax
2758	read (ilesfile,*) height(k),thlprof(k),qtprof (k), &
2759	& uprof (k),vprof (k),e12prof(k)
2760	enddo
2761	close(ilesfile)
2762
2763	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2764	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2765	read (ilesfile,*) kmax2
2766	if (kmax .ne. kmax2) then
2767	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2768	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2769	stop 'lecture profiles'
2770	endif
2771	do k=1,kmax
2772	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &
2773	& dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2774	end do
2775	close(ilesfile)
2776
2777	open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
2778	if (ierr /= 0) then
2779	print*,'WARNING : trac.inp does not exist'
2780	else
2781	read (ilesfile,*) kmax2,nt1,nt2
2782	if (nt2>ntrac) then
2783	stop'Augmenter le nombre de traceurs dans traceur.def'
2784	endif
2785	if (kmax .ne. kmax2) then
2786	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2787	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2788	stop 'lecture profiles'
2789	endif
2790	do k=1,kmax
2791	read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
2792	end do
2793	close(ilesfile)
2794	endif
2795
2796	return
2797	end
2798	!======================================================================
2799	subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &
2800	& thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
2801	!======================================================================
2802	implicit none
2803
2804	integer nlev_max,kmax
2805	logical :: llesread = .true.
2806
2807	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2808	real thlprof(nlev_max)
2809	real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2810	real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
2811
2812	integer, parameter :: ilesfile=1
2813	integer :: k,ierr
2814
2815	if(.not.(llesread)) return
2816
2817	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2818	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2819	read (ilesfile,*) kmax
2820	do k=1,kmax
2821	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
2822	& qprof (k),uprof(k), vprof(k), wprof(k), &
2823	& omega (k),o3mmr(k)
2824	enddo
2825	close(ilesfile)
2826
2827	return
2828	end
2829
2830	!======================================================================
2831	subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &
2832	& thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
2833	!======================================================================
2834	implicit none
2835
2836	integer nlev_max,kmax
2837	logical :: llesread = .true.
2838
2839	real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &
2840	& thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &
2841	& qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &
2842	& wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
2843
2844	integer, parameter :: ilesfile=1
2845	integer :: ierr,k
2846
2847	if(.not.(llesread)) return
2848
2849	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2850	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2851	read (ilesfile,*) kmax
2852	do k=1,kmax
2853	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
2854	& qvprof (k),qlprof (k),qtprof (k), &
2855	& uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
2856	enddo
2857	close(ilesfile)
2858
2859	return
2860	end
2861
2862
2863
2864	!======================================================================
2865	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &
2866	& vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2867	!======================================================================
2868	implicit none
2869
2870	integer nlev_max,kmax
2871	logical :: llesread = .true.
2872
2873	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2874	real thetaprof(nlev_max),rvprof(nlev_max)
2875	real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2876	real aprof(nlev_max+1),bprof(nlev_max+1)
2877
2878	integer, parameter :: ilesfile=1
2879	integer, parameter :: ifile=2
2880	integer :: ierr,jtot,k
2881
2882	if(.not.(llesread)) return
2883
2884	! Read profiles at full levels
2885	IF(nlev_max.EQ.19) THEN
2886	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2887	print *,'On ouvre prof.inp.19'
2888	ELSE
2889	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2890	print *,'On ouvre prof.inp.40'
2891	ENDIF
2892	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2893	read (ilesfile,*) kmax
2894	do k=1,kmax
2895	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &
2896	& thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2897	enddo
2898	close(ilesfile)
2899
2900	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
2901	IF(nlev_max.EQ.19) THEN
2902	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2903	print *,'On ouvre proh.inp.19'
2904	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2905	ELSE
2906	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2907	print *,'On ouvre proh.inp.40'
2908	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2909	ENDIF
2910	read (ifile,*) kmax
2911	do k=1,kmax
2912	read (ifile,*) jtot,aprof(k),bprof(k)
2913	enddo
2914	close(ifile)
2915
2916	return
2917	end
2918	!=====================================================================
2919	subroutine read_amma(fich_amma,nlevel,ntime &
2920	& ,zz,pp,temp,qv,u,v,dw &
2921	& ,dt,dq,sens,flat)
2922
2923	!program reading forcings of the AMMA case study
2924
2925
2926	implicit none
2927
2928	#include "netcdf.inc"
2929
2930	integer ntime,nlevel
2931	character*80 :: fich_amma
2932	real*8 zz(nlevel)
2933
2934	real*8 temp(nlevel),pp(nlevel)
2935	real*8 qv(nlevel),u(nlevel)
2936	real*8 v(nlevel)
2937	real*8 dw(nlevel,ntime)
2938	real*8 dt(nlevel,ntime)
2939	real*8 dq(nlevel,ntime)
2940	real*8 flat(ntime),sens(ntime)
2941
2942	integer nid, ierr
2943	integer nbvar3d
2944	parameter(nbvar3d=30)
2945	integer var3didin(nbvar3d)
2946
2947	ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
2948	if (ierr.NE.NF_NOERR) then
2949	write(,) 'ERROR: Pb opening forcings nc file '
2950	write(,) NF_STRERROR(ierr)
2951	stop ""
2952	endif
2953
2954
2955	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
2956	if(ierr/=NF_NOERR) then
2957	write(,) NF_STRERROR(ierr)
2958	stop 'lev'
2959	endif
2960
2961
2962	ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
2963	if(ierr/=NF_NOERR) then
2964	write(,) NF_STRERROR(ierr)
2965	stop 'temp'
2966	endif
2967
2968	ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
2969	if(ierr/=NF_NOERR) then
2970	write(,) NF_STRERROR(ierr)
2971	stop 'qv'
2972	endif
2973
2974	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
2975	if(ierr/=NF_NOERR) then
2976	write(,) NF_STRERROR(ierr)
2977	stop 'u'
2978	endif
2979
2980	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
2981	if(ierr/=NF_NOERR) then
2982	write(,) NF_STRERROR(ierr)
2983	stop 'v'
2984	endif
2985
2986	ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
2987	if(ierr/=NF_NOERR) then
2988	write(,) NF_STRERROR(ierr)
2989	stop 'dw'
2990	endif
2991
2992	ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
2993	if(ierr/=NF_NOERR) then
2994	write(,) NF_STRERROR(ierr)
2995	stop 'dt'
2996	endif
2997
2998	ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
2999	if(ierr/=NF_NOERR) then
3000	write(,) NF_STRERROR(ierr)
3001	stop 'dq'
3002	endif
3003
3004	ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
3005	if(ierr/=NF_NOERR) then
3006	write(,) NF_STRERROR(ierr)
3007	stop 'sens'
3008	endif
3009
3010	ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
3011	if(ierr/=NF_NOERR) then
3012	write(,) NF_STRERROR(ierr)
3013	stop 'flat'
3014	endif
3015
3016	ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
3017	if(ierr/=NF_NOERR) then
3018	write(,) NF_STRERROR(ierr)
3019	stop 'pp'
3020	endif
3021
3022	!dimensions lecture
3023	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3024
3025	#ifdef NC_DOUBLE
3026	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3027	#else
3028	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3029	#endif
3030	if(ierr/=NF_NOERR) then
3031	write(,) NF_STRERROR(ierr)
3032	stop "getvarup"
3033	endif
3034	! write(,)'lecture z ok',zz
3035
3036	#ifdef NC_DOUBLE
3037	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
3038	#else
3039	ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
3040	#endif
3041	if(ierr/=NF_NOERR) then
3042	write(,) NF_STRERROR(ierr)
3043	stop "getvarup"
3044	endif
3045	! write(,)'lecture th ok',temp
3046
3047	#ifdef NC_DOUBLE
3048	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
3049	#else
3050	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
3051	#endif
3052	if(ierr/=NF_NOERR) then
3053	write(,) NF_STRERROR(ierr)
3054	stop "getvarup"
3055	endif
3056	! write(,)'lecture qv ok',qv
3057
3058	#ifdef NC_DOUBLE
3059	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3060	#else
3061	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3062	#endif
3063	if(ierr/=NF_NOERR) then
3064	write(,) NF_STRERROR(ierr)
3065	stop "getvarup"
3066	endif
3067	! write(,)'lecture u ok',u
3068
3069	#ifdef NC_DOUBLE
3070	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3071	#else
3072	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3073	#endif
3074	if(ierr/=NF_NOERR) then
3075	write(,) NF_STRERROR(ierr)
3076	stop "getvarup"
3077	endif
3078	! write(,)'lecture v ok',v
3079
3080	#ifdef NC_DOUBLE
3081	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
3082	#else
3083	ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
3084	#endif
3085	if(ierr/=NF_NOERR) then
3086	write(,) NF_STRERROR(ierr)
3087	stop "getvarup"
3088	endif
3089	! write(,)'lecture w ok',dw
3090
3091	#ifdef NC_DOUBLE
3092	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
3093	#else
3094	ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
3095	#endif
3096	if(ierr/=NF_NOERR) then
3097	write(,) NF_STRERROR(ierr)
3098	stop "getvarup"
3099	endif
3100	! write(,)'lecture dt ok',dt
3101
3102	#ifdef NC_DOUBLE
3103	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
3104	#else
3105	ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
3106	#endif
3107	if(ierr/=NF_NOERR) then
3108	write(,) NF_STRERROR(ierr)
3109	stop "getvarup"
3110	endif
3111	! write(,)'lecture dq ok',dq
3112
3113	#ifdef NC_DOUBLE
3114	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
3115	#else
3116	ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
3117	#endif
3118	if(ierr/=NF_NOERR) then
3119	write(,) NF_STRERROR(ierr)
3120	stop "getvarup"
3121	endif
3122	! write(,)'lecture sens ok',sens
3123
3124	#ifdef NC_DOUBLE
3125	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
3126	#else
3127	ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
3128	#endif
3129	if(ierr/=NF_NOERR) then
3130	write(,) NF_STRERROR(ierr)
3131	stop "getvarup"
3132	endif
3133	! write(,)'lecture flat ok',flat
3134
3135	#ifdef NC_DOUBLE
3136	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
3137	#else
3138	ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
3139	#endif
3140	if(ierr/=NF_NOERR) then
3141	write(,) NF_STRERROR(ierr)
3142	stop "getvarup"
3143	endif
3144	! write(,)'lecture pp ok',pp
3145
3146	return
3147	end subroutine read_amma
3148	!======================================================================
3149	SUBROUTINE interp_amma_time(day,day1,annee_ref &
3150	& ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma &
3151	& ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma &
3152	& ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
3153	implicit none
3154
3155	!---------------------------------------------------------------------------------------
3156	! Time interpolation of a 2D field to the timestep corresponding to day
3157	!
3158	! day: current julian day (e.g. 717538.2)
3159	! day1: first day of the simulation
3160	! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
3161	! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
3162	!---------------------------------------------------------------------------------------
3163
3164	#include "compar1d.h"
3165
3166	! inputs:
3167	integer annee_ref
3168	integer nt_amma,nlev_amma
3169	integer year_ini_amma
3170	real day, day1,day_ini_amma,dt_amma
3171	real vitw_amma(nlev_amma,nt_amma)
3172	real ht_amma(nlev_amma,nt_amma)
3173	real hq_amma(nlev_amma,nt_amma)
3174	real lat_amma(nt_amma)
3175	real sens_amma(nt_amma)
3176	! outputs:
3177	real vitw_prof(nlev_amma)
3178	real ht_prof(nlev_amma)
3179	real hq_prof(nlev_amma)
3180	real lat_prof,sens_prof
3181	! local:
3182	integer it_amma1, it_amma2,k
3183	real timeit,time_amma1,time_amma2,frac
3184
3185
3186	if (forcing_type.eq.6) then
3187	! Check that initial day of the simulation consistent with AMMA case:
3188	if (annee_ref.ne.2006) then
3189	print*,'Pour AMMA, annee_ref doit etre 2006'
3190	print*,'Changer annee_ref dans run.def'
3191	stop
3192	endif
3193	if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
3194	print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
3195	print*,'Changer dayref dans run.def'
3196	stop
3197	endif
3198	if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
3199	print*,'AMMA a fini le 11 juillet'
3200	print*,'Changer dayref ou nday dans run.def'
3201	stop
3202	endif
3203	endif
3204
3205	! Determine timestep relative to the 1st day of AMMA:
3206	! timeit=(day-day1)*86400.
3207	! if (annee_ref.eq.1992) then
3208	! timeit=(day-day_ini_toga)*86400.
3209	! else
3210	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3211	! endif
3212	timeit=(day-day_ini_amma)*86400
3213
3214	! Determine the closest observation times:
3215	! it_amma1=INT(timeit/dt_amma)+1
3216	! it_amma2=it_amma1 + 1
3217	! time_amma1=(it_amma1-1)*dt_amma
3218	! time_amma2=(it_amma2-1)*dt_amma
3219
3220	it_amma1=INT(timeit/dt_amma)+1
3221	IF (it_amma1 .EQ. nt_amma) THEN
3222	it_amma2=it_amma1
3223	ELSE
3224	it_amma2=it_amma1 + 1
3225	ENDIF
3226	time_amma1=(it_amma1-1)*dt_amma
3227	time_amma2=(it_amma2-1)*dt_amma
3228
3229	if (it_amma1 .gt. nt_amma) then
3230	write(,) 'PB-stop: day, it_amma1, it_amma2, timeit: ' &
3231	& ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
3232	stop
3233	endif
3234
3235	! time interpolation:
3236	frac=(time_amma2-timeit)/(time_amma2-time_amma1)
3237	frac=max(frac,0.0)
3238
3239	lat_prof = lat_amma(it_amma2) &
3240	& -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
3241	sens_prof = sens_amma(it_amma2) &
3242	& -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
3243
3244	do k=1,nlev_amma
3245	vitw_prof(k) = vitw_amma(k,it_amma2) &
3246	& -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
3247	ht_prof(k) = ht_amma(k,it_amma2) &
3248	& -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
3249	hq_prof(k) = hq_amma(k,it_amma2) &
3250	& -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
3251	enddo
3252
3253	return
3254	END
3255
3256	!=====================================================================
3257	subroutine read_fire(fich_fire,nlevel,ntime &
3258	& ,zz,thl,qt,u,v,tke &
3259	& ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3260
3261	!program reading forcings of the FIRE case study
3262
3263
3264	implicit none
3265
3266	#include "netcdf.inc"
3267
3268	integer ntime,nlevel
3269	character*80 :: fich_fire
3270	real*8 zz(nlevel)
3271
3272	real*8 thl(nlevel)
3273	real*8 qt(nlevel),u(nlevel)
3274	real*8 v(nlevel),tke(nlevel)
3275	real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3276	real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3277	real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3278
3279	integer nid, ierr
3280	integer nbvar3d
3281	parameter(nbvar3d=30)
3282	integer var3didin(nbvar3d)
3283
3284	ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3285	if (ierr.NE.NF_NOERR) then
3286	write(,) 'ERROR: Pb opening forcings nc file '
3287	write(,) NF_STRERROR(ierr)
3288	stop ""
3289	endif
3290
3291
3292	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
3293	if(ierr/=NF_NOERR) then
3294	write(,) NF_STRERROR(ierr)
3295	stop 'lev'
3296	endif
3297
3298
3299	ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3300	if(ierr/=NF_NOERR) then
3301	write(,) NF_STRERROR(ierr)
3302	stop 'temp'
3303	endif
3304
3305	ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3306	if(ierr/=NF_NOERR) then
3307	write(,) NF_STRERROR(ierr)
3308	stop 'qv'
3309	endif
3310
3311	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3312	if(ierr/=NF_NOERR) then
3313	write(,) NF_STRERROR(ierr)
3314	stop 'u'
3315	endif
3316
3317	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3318	if(ierr/=NF_NOERR) then
3319	write(,) NF_STRERROR(ierr)
3320	stop 'v'
3321	endif
3322
3323	ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3324	if(ierr/=NF_NOERR) then
3325	write(,) NF_STRERROR(ierr)
3326	stop 'tke'
3327	endif
3328
3329	ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3330	if(ierr/=NF_NOERR) then
3331	write(,) NF_STRERROR(ierr)
3332	stop 'ug'
3333	endif
3334
3335	ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3336	if(ierr/=NF_NOERR) then
3337	write(,) NF_STRERROR(ierr)
3338	stop 'vg'
3339	endif
3340
3341	ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3342	if(ierr/=NF_NOERR) then
3343	write(,) NF_STRERROR(ierr)
3344	stop 'wls'
3345	endif
3346
3347	ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3348	if(ierr/=NF_NOERR) then
3349	write(,) NF_STRERROR(ierr)
3350	stop 'dqtdx'
3351	endif
3352
3353	ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3354	if(ierr/=NF_NOERR) then
3355	write(,) NF_STRERROR(ierr)
3356	stop 'dqtdy'
3357	endif
3358
3359	ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3360	if(ierr/=NF_NOERR) then
3361	write(,) NF_STRERROR(ierr)
3362	stop 'dqtdt'
3363	endif
3364
3365	ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3366	if(ierr/=NF_NOERR) then
3367	write(,) NF_STRERROR(ierr)
3368	stop 'thl_rad'
3369	endif
3370	!dimensions lecture
3371	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3372
3373	#ifdef NC_DOUBLE
3374	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3375	#else
3376	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3377	#endif
3378	if(ierr/=NF_NOERR) then
3379	write(,) NF_STRERROR(ierr)
3380	stop "getvarup"
3381	endif
3382	! write(,)'lecture z ok',zz
3383
3384	#ifdef NC_DOUBLE
3385	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3386	#else
3387	ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3388	#endif
3389	if(ierr/=NF_NOERR) then
3390	write(,) NF_STRERROR(ierr)
3391	stop "getvarup"
3392	endif
3393	! write(,)'lecture thl ok',thl
3394
3395	#ifdef NC_DOUBLE
3396	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3397	#else
3398	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3399	#endif
3400	if(ierr/=NF_NOERR) then
3401	write(,) NF_STRERROR(ierr)
3402	stop "getvarup"
3403	endif
3404	! write(,)'lecture qt ok',qt
3405
3406	#ifdef NC_DOUBLE
3407	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3408	#else
3409	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3410	#endif
3411	if(ierr/=NF_NOERR) then
3412	write(,) NF_STRERROR(ierr)
3413	stop "getvarup"
3414	endif
3415	! write(,)'lecture u ok',u
3416
3417	#ifdef NC_DOUBLE
3418	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3419	#else
3420	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3421	#endif
3422	if(ierr/=NF_NOERR) then
3423	write(,) NF_STRERROR(ierr)
3424	stop "getvarup"
3425	endif
3426	! write(,)'lecture v ok',v
3427
3428	#ifdef NC_DOUBLE
3429	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3430	#else
3431	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3432	#endif
3433	if(ierr/=NF_NOERR) then
3434	write(,) NF_STRERROR(ierr)
3435	stop "getvarup"
3436	endif
3437	! write(,)'lecture tke ok',tke
3438
3439	#ifdef NC_DOUBLE
3440	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3441	#else
3442	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3443	#endif
3444	if(ierr/=NF_NOERR) then
3445	write(,) NF_STRERROR(ierr)
3446	stop "getvarup"
3447	endif
3448	! write(,)'lecture ug ok',ug
3449
3450	#ifdef NC_DOUBLE
3451	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3452	#else
3453	ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3454	#endif
3455	if(ierr/=NF_NOERR) then
3456	write(,) NF_STRERROR(ierr)
3457	stop "getvarup"
3458	endif
3459	! write(,)'lecture vg ok',vg
3460
3461	#ifdef NC_DOUBLE
3462	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3463	#else
3464	ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3465	#endif
3466	if(ierr/=NF_NOERR) then
3467	write(,) NF_STRERROR(ierr)
3468	stop "getvarup"
3469	endif
3470	! write(,)'lecture wls ok',wls
3471
3472	#ifdef NC_DOUBLE
3473	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3474	#else
3475	ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3476	#endif
3477	if(ierr/=NF_NOERR) then
3478	write(,) NF_STRERROR(ierr)
3479	stop "getvarup"
3480	endif
3481	! write(,)'lecture dqtdx ok',dqtdx
3482
3483	#ifdef NC_DOUBLE
3484	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3485	#else
3486	ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3487	#endif
3488	if(ierr/=NF_NOERR) then
3489	write(,) NF_STRERROR(ierr)
3490	stop "getvarup"
3491	endif
3492	! write(,)'lecture dqtdy ok',dqtdy
3493
3494	#ifdef NC_DOUBLE
3495	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3496	#else
3497	ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3498	#endif
3499	if(ierr/=NF_NOERR) then
3500	write(,) NF_STRERROR(ierr)
3501	stop "getvarup"
3502	endif
3503	! write(,)'lecture dqtdt ok',dqtdt
3504
3505	#ifdef NC_DOUBLE
3506	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3507	#else
3508	ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3509	#endif
3510	if(ierr/=NF_NOERR) then
3511	write(,) NF_STRERROR(ierr)
3512	stop "getvarup"
3513	endif
3514	! write(,)'lecture thl_rad ok',thl_rad
3515
3516	return
3517	end subroutine read_fire

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

Download in other formats: