Context Navigation

source: LMDZ5/trunk/libf/phy1d/1DUTILS.h @ 1975

Last change on this file since 1975 was 1975, checked in by fhourdin, 10 years ago
On enleve wrgradsfi de 1DUTILS.h. La routine est deja dans phylmd. Removing wrgradsfi from 1DUTILS.h (already in phylmd)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 109.8 KB

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

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

Download in other formats: