Context Navigation

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

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

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

Download in other formats: