Context Navigation

source: LMDZ5/branches/testing/libf/phy1d/1DUTILS.h @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 108.1 KB

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

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

Download in other formats: