Context Navigation

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

Last change on this file since 2019 was 2019, checked in by fhourdin, 10 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: