Context Navigation

source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/1DUTILS.h

Last change on this file was 2019, checked in by fhourdin, 11 years ago
Passage au format libre pour inclure les ficheirs du 1D dans lmdz1d.F90
File size: 111.4 KB

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

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

Download in other formats: