Context Navigation

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

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

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

Download in other formats: