Context Navigation

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

Last change on this file since 1953 was 1953, checked in by fhourdin, 11 years ago
Correction in for readprofiles_fire
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	if (ierr /= 0) then
3020	print*,'WARNING : trac.inp does not exist'
3021	else
3022	read (ilesfile,*) kmax2,nt1,nt2
3023	if (nt2>ntrac) then
3024	stop'Augmenter le nombre de traceurs dans traceur.def'
3025	endif
3026	if (kmax .ne. kmax2) then
3027	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3028	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3029	stop 'lecture profiles'
3030	endif
3031	do k=1,kmax
3032	read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
3033	end do
3034	close(ilesfile)
3035	endif
3036
3037	return
3038	end
3039	!======================================================================
3040	subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,
3041	. thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
3042	!======================================================================
3043	implicit none
3044
3045	integer nlev_max,kmax,kmax2
3046	logical :: llesread = .true.
3047
3048	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
3049	. thlprof(nlev_max),
3050	. qprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
3051	. wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
3052
3053	integer, parameter :: ilesfile=1
3054	integer :: ierr,irad,imax,jtot,k
3055	logical :: lmoist,lcoriol,ltimedep
3056	real :: xsize,ysize
3057	real :: ustin,wsvsurf,timerad
3058	character(80) :: chmess
3059
3060	if(.not.(llesread)) return
3061
3062	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3063	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3064	read (ilesfile,*) kmax
3065	do k=1,kmax
3066	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k),
3067	. qprof (k),uprof(k), vprof(k), wprof(k),
3068	. omega (k),o3mmr(k)
3069	enddo
3070	close(ilesfile)
3071
3072	return
3073	end
3074
3075	!======================================================================
3076	subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,
3077	. thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
3078	!======================================================================
3079	implicit none
3080
3081	integer nlev_max,kmax,kmax2
3082	logical :: llesread = .true.
3083
3084	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
3085	. thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),
3086	. qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
3087	. wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
3088
3089	integer, parameter :: ilesfile=1
3090	integer :: ierr,irad,imax,jtot,k
3091	logical :: lmoist,lcoriol,ltimedep
3092	real :: xsize,ysize
3093	real :: ustin,wsvsurf,timerad
3094	character(80) :: chmess
3095
3096	if(.not.(llesread)) return
3097
3098	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3099	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3100	read (ilesfile,*) kmax
3101	do k=1,kmax
3102	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k),
3103	. qvprof (k),qlprof (k),qtprof (k),
3104	. uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
3105	enddo
3106	close(ilesfile)
3107
3108	return
3109	end
3110
3111
3112
3113	!======================================================================
3114	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
3115	. vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
3116	!======================================================================
3117	implicit none
3118
3119	integer nlev_max,kmax,kmax2
3120	logical :: llesread = .true.
3121
3122	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
3123	. thetaprof(nlev_max),rvprof(nlev_max),
3124	. qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
3125	. aprof(nlev_max+1),bprof(nlev_max+1)
3126
3127	integer, parameter :: ilesfile=1
3128	integer, parameter :: ifile=2
3129	integer :: ierr,irad,imax,jtot,k
3130	logical :: lmoist,lcoriol,ltimedep
3131	real :: xsize,ysize
3132	real :: ustin,wsvsurf,timerad
3133	character(80) :: chmess
3134
3135	if(.not.(llesread)) return
3136
3137	! Read profiles at full levels
3138	IF(nlev_max.EQ.19) THEN
3139	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
3140	print *,'On ouvre prof.inp.19'
3141	ELSE
3142	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
3143	print *,'On ouvre prof.inp.40'
3144	ENDIF
3145	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3146	read (ilesfile,*) kmax
3147	do k=1,kmax
3148	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k),
3149	. thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
3150	enddo
3151	close(ilesfile)
3152
3153	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
3154	IF(nlev_max.EQ.19) THEN
3155	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
3156	print *,'On ouvre proh.inp.19'
3157	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
3158	ELSE
3159	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
3160	print *,'On ouvre proh.inp.40'
3161	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
3162	ENDIF
3163	read (ifile,*) kmax
3164	do k=1,kmax
3165	read (ifile,*) jtot,aprof(k),bprof(k)
3166	enddo
3167	close(ifile)
3168
3169	return
3170	end
3171	!=====================================================================
3172	subroutine read_amma(fich_amma,nlevel,ntime
3173	: ,zz,pp,temp,qv,u,v,dw
3174	: ,dt,dq,sens,flat)
3175
3176	!program reading forcings of the AMMA case study
3177
3178
3179	implicit none
3180
3181	#include "netcdf.inc"
3182
3183	integer ntime,nlevel
3184	integer l,k
3185	character*80 :: fich_amma
3186	real*8 time(ntime)
3187	real*8 zz(nlevel)
3188
3189	real*8 temp(nlevel),pp(nlevel)
3190	real*8 qv(nlevel),u(nlevel)
3191	real*8 v(nlevel)
3192	real*8 dw(nlevel,ntime)
3193	real*8 dt(nlevel,ntime)
3194	real*8 dq(nlevel,ntime)
3195	real*8 flat(ntime),sens(ntime)
3196
3197	integer nid, ierr
3198	integer nbvar3d
3199	parameter(nbvar3d=30)
3200	integer var3didin(nbvar3d)
3201
3202	ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
3203	if (ierr.NE.NF_NOERR) then
3204	write(,) 'ERROR: Pb opening forcings nc file '
3205	write(,) NF_STRERROR(ierr)
3206	stop ""
3207	endif
3208
3209
3210	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
3211	if(ierr/=NF_NOERR) then
3212	write(,) NF_STRERROR(ierr)
3213	stop 'lev'
3214	endif
3215
3216
3217	ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
3218	if(ierr/=NF_NOERR) then
3219	write(,) NF_STRERROR(ierr)
3220	stop 'temp'
3221	endif
3222
3223	ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
3224	if(ierr/=NF_NOERR) then
3225	write(,) NF_STRERROR(ierr)
3226	stop 'qv'
3227	endif
3228
3229	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3230	if(ierr/=NF_NOERR) then
3231	write(,) NF_STRERROR(ierr)
3232	stop 'u'
3233	endif
3234
3235	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3236	if(ierr/=NF_NOERR) then
3237	write(,) NF_STRERROR(ierr)
3238	stop 'v'
3239	endif
3240
3241	ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
3242	if(ierr/=NF_NOERR) then
3243	write(,) NF_STRERROR(ierr)
3244	stop 'dw'
3245	endif
3246
3247	ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
3248	if(ierr/=NF_NOERR) then
3249	write(,) NF_STRERROR(ierr)
3250	stop 'dt'
3251	endif
3252
3253	ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
3254	if(ierr/=NF_NOERR) then
3255	write(,) NF_STRERROR(ierr)
3256	stop 'dq'
3257	endif
3258
3259	ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
3260	if(ierr/=NF_NOERR) then
3261	write(,) NF_STRERROR(ierr)
3262	stop 'sens'
3263	endif
3264
3265	ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
3266	if(ierr/=NF_NOERR) then
3267	write(,) NF_STRERROR(ierr)
3268	stop 'flat'
3269	endif
3270
3271	ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
3272	if(ierr/=NF_NOERR) then
3273	write(,) NF_STRERROR(ierr)
3274	stop 'pp'
3275	endif
3276
3277	!dimensions lecture
3278	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3279
3280	#ifdef NC_DOUBLE
3281	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3282	#else
3283	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3284	#endif
3285	if(ierr/=NF_NOERR) then
3286	write(,) NF_STRERROR(ierr)
3287	stop "getvarup"
3288	endif
3289	! write(,)'lecture z ok',zz
3290
3291	#ifdef NC_DOUBLE
3292	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
3293	#else
3294	ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
3295	#endif
3296	if(ierr/=NF_NOERR) then
3297	write(,) NF_STRERROR(ierr)
3298	stop "getvarup"
3299	endif
3300	! write(,)'lecture th ok',temp
3301
3302	#ifdef NC_DOUBLE
3303	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
3304	#else
3305	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
3306	#endif
3307	if(ierr/=NF_NOERR) then
3308	write(,) NF_STRERROR(ierr)
3309	stop "getvarup"
3310	endif
3311	! write(,)'lecture qv ok',qv
3312
3313	#ifdef NC_DOUBLE
3314	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3315	#else
3316	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3317	#endif
3318	if(ierr/=NF_NOERR) then
3319	write(,) NF_STRERROR(ierr)
3320	stop "getvarup"
3321	endif
3322	! write(,)'lecture u ok',u
3323
3324	#ifdef NC_DOUBLE
3325	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3326	#else
3327	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3328	#endif
3329	if(ierr/=NF_NOERR) then
3330	write(,) NF_STRERROR(ierr)
3331	stop "getvarup"
3332	endif
3333	! write(,)'lecture v ok',v
3334
3335	#ifdef NC_DOUBLE
3336	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
3337	#else
3338	ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
3339	#endif
3340	if(ierr/=NF_NOERR) then
3341	write(,) NF_STRERROR(ierr)
3342	stop "getvarup"
3343	endif
3344	! write(,)'lecture w ok',dw
3345
3346	#ifdef NC_DOUBLE
3347	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
3348	#else
3349	ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
3350	#endif
3351	if(ierr/=NF_NOERR) then
3352	write(,) NF_STRERROR(ierr)
3353	stop "getvarup"
3354	endif
3355	! write(,)'lecture dt ok',dt
3356
3357	#ifdef NC_DOUBLE
3358	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
3359	#else
3360	ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
3361	#endif
3362	if(ierr/=NF_NOERR) then
3363	write(,) NF_STRERROR(ierr)
3364	stop "getvarup"
3365	endif
3366	! write(,)'lecture dq ok',dq
3367
3368	#ifdef NC_DOUBLE
3369	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
3370	#else
3371	ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
3372	#endif
3373	if(ierr/=NF_NOERR) then
3374	write(,) NF_STRERROR(ierr)
3375	stop "getvarup"
3376	endif
3377	! write(,)'lecture sens ok',sens
3378
3379	#ifdef NC_DOUBLE
3380	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
3381	#else
3382	ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
3383	#endif
3384	if(ierr/=NF_NOERR) then
3385	write(,) NF_STRERROR(ierr)
3386	stop "getvarup"
3387	endif
3388	! write(,)'lecture flat ok',flat
3389
3390	#ifdef NC_DOUBLE
3391	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
3392	#else
3393	ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
3394	#endif
3395	if(ierr/=NF_NOERR) then
3396	write(,) NF_STRERROR(ierr)
3397	stop "getvarup"
3398	endif
3399	! write(,)'lecture pp ok',pp
3400
3401	return
3402	end subroutine read_amma
3403	!======================================================================
3404	SUBROUTINE interp_amma_time(day,day1,annee_ref
3405	i ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma
3406	i ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
3407	o ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
3408	implicit none
3409
3410	!---------------------------------------------------------------------------------------
3411	! Time interpolation of a 2D field to the timestep corresponding to day
3412	!
3413	! day: current julian day (e.g. 717538.2)
3414	! day1: first day of the simulation
3415	! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
3416	! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
3417	!---------------------------------------------------------------------------------------
3418
3419	#include "compar1d.h"
3420
3421	! inputs:
3422	integer annee_ref
3423	integer nt_amma,nlev_amma
3424	integer year_ini_amma
3425	real day, day1,day_ini_amma,dt_amma
3426	real vitw_amma(nlev_amma,nt_amma)
3427	real ht_amma(nlev_amma,nt_amma)
3428	real hq_amma(nlev_amma,nt_amma)
3429	real lat_amma(nt_amma)
3430	real sens_amma(nt_amma)
3431	! outputs:
3432	real vitw_prof(nlev_amma)
3433	real ht_prof(nlev_amma)
3434	real hq_prof(nlev_amma)
3435	real lat_prof,sens_prof
3436	! local:
3437	integer it_amma1, it_amma2,k
3438	real timeit,time_amma1,time_amma2,frac
3439
3440
3441	if (forcing_type.eq.6) then
3442	! Check that initial day of the simulation consistent with AMMA case:
3443	if (annee_ref.ne.2006) then
3444	print*,'Pour AMMA, annee_ref doit etre 2006'
3445	print*,'Changer annee_ref dans run.def'
3446	stop
3447	endif
3448	if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
3449	print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
3450	print*,'Changer dayref dans run.def'
3451	stop
3452	endif
3453	if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
3454	print*,'AMMA a fini le 11 juillet'
3455	print*,'Changer dayref ou nday dans run.def'
3456	stop
3457	endif
3458	endif
3459
3460	! Determine timestep relative to the 1st day of AMMA:
3461	! timeit=(day-day1)*86400.
3462	! if (annee_ref.eq.1992) then
3463	! timeit=(day-day_ini_toga)*86400.
3464	! else
3465	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3466	! endif
3467	timeit=(day-day_ini_amma)*86400
3468
3469	! Determine the closest observation times:
3470	! it_amma1=INT(timeit/dt_amma)+1
3471	! it_amma2=it_amma1 + 1
3472	! time_amma1=(it_amma1-1)*dt_amma
3473	! time_amma2=(it_amma2-1)*dt_amma
3474
3475	it_amma1=INT(timeit/dt_amma)+1
3476	IF (it_amma1 .EQ. nt_amma) THEN
3477	it_amma2=it_amma1
3478	ELSE
3479	it_amma2=it_amma1 + 1
3480	ENDIF
3481	time_amma1=(it_amma1-1)*dt_amma
3482	time_amma2=(it_amma2-1)*dt_amma
3483
3484	if (it_amma1 .gt. nt_amma) then
3485	write(,) 'PB-stop: day, it_amma1, it_amma2, timeit: '
3486	: ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
3487	stop
3488	endif
3489
3490	! time interpolation:
3491	frac=(time_amma2-timeit)/(time_amma2-time_amma1)
3492	frac=max(frac,0.0)
3493
3494	lat_prof = lat_amma(it_amma2)
3495	: -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
3496	sens_prof = sens_amma(it_amma2)
3497	: -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
3498
3499	do k=1,nlev_amma
3500	vitw_prof(k) = vitw_amma(k,it_amma2)
3501	: -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
3502	ht_prof(k) = ht_amma(k,it_amma2)
3503	: -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
3504	hq_prof(k) = hq_amma(k,it_amma2)
3505	: -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
3506	enddo
3507
3508	return
3509	END
3510
3511	!=====================================================================
3512	subroutine read_fire(fich_fire,nlevel,ntime
3513	: ,zz,thl,qt,u,v,tke
3514	: ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3515
3516	!program reading forcings of the FIRE case study
3517
3518
3519	implicit none
3520
3521	#include "netcdf.inc"
3522
3523	integer ntime,nlevel
3524	integer l,k
3525	character*80 :: fich_fire
3526	real*8 time(ntime)
3527	real*8 zz(nlevel)
3528
3529	real*8 thl(nlevel)
3530	real*8 qt(nlevel),u(nlevel)
3531	real*8 v(nlevel),tke(nlevel)
3532	real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3533	real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3534	real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3535
3536	integer nid, ierr
3537	integer nbvar3d
3538	parameter(nbvar3d=30)
3539	integer var3didin(nbvar3d)
3540
3541	ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3542	if (ierr.NE.NF_NOERR) then
3543	write(,) 'ERROR: Pb opening forcings nc file '
3544	write(,) NF_STRERROR(ierr)
3545	stop ""
3546	endif
3547
3548
3549	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
3550	if(ierr/=NF_NOERR) then
3551	write(,) NF_STRERROR(ierr)
3552	stop 'lev'
3553	endif
3554
3555
3556	ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3557	if(ierr/=NF_NOERR) then
3558	write(,) NF_STRERROR(ierr)
3559	stop 'temp'
3560	endif
3561
3562	ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3563	if(ierr/=NF_NOERR) then
3564	write(,) NF_STRERROR(ierr)
3565	stop 'qv'
3566	endif
3567
3568	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3569	if(ierr/=NF_NOERR) then
3570	write(,) NF_STRERROR(ierr)
3571	stop 'u'
3572	endif
3573
3574	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3575	if(ierr/=NF_NOERR) then
3576	write(,) NF_STRERROR(ierr)
3577	stop 'v'
3578	endif
3579
3580	ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3581	if(ierr/=NF_NOERR) then
3582	write(,) NF_STRERROR(ierr)
3583	stop 'tke'
3584	endif
3585
3586	ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3587	if(ierr/=NF_NOERR) then
3588	write(,) NF_STRERROR(ierr)
3589	stop 'ug'
3590	endif
3591
3592	ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3593	if(ierr/=NF_NOERR) then
3594	write(,) NF_STRERROR(ierr)
3595	stop 'vg'
3596	endif
3597
3598	ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3599	if(ierr/=NF_NOERR) then
3600	write(,) NF_STRERROR(ierr)
3601	stop 'wls'
3602	endif
3603
3604	ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3605	if(ierr/=NF_NOERR) then
3606	write(,) NF_STRERROR(ierr)
3607	stop 'dqtdx'
3608	endif
3609
3610	ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3611	if(ierr/=NF_NOERR) then
3612	write(,) NF_STRERROR(ierr)
3613	stop 'dqtdy'
3614	endif
3615
3616	ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3617	if(ierr/=NF_NOERR) then
3618	write(,) NF_STRERROR(ierr)
3619	stop 'dqtdt'
3620	endif
3621
3622	ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3623	if(ierr/=NF_NOERR) then
3624	write(,) NF_STRERROR(ierr)
3625	stop 'thl_rad'
3626	endif
3627	!dimensions lecture
3628	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3629
3630	#ifdef NC_DOUBLE
3631	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3632	#else
3633	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3634	#endif
3635	if(ierr/=NF_NOERR) then
3636	write(,) NF_STRERROR(ierr)
3637	stop "getvarup"
3638	endif
3639	! write(,)'lecture z ok',zz
3640
3641	#ifdef NC_DOUBLE
3642	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3643	#else
3644	ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3645	#endif
3646	if(ierr/=NF_NOERR) then
3647	write(,) NF_STRERROR(ierr)
3648	stop "getvarup"
3649	endif
3650	! write(,)'lecture thl ok',thl
3651
3652	#ifdef NC_DOUBLE
3653	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3654	#else
3655	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3656	#endif
3657	if(ierr/=NF_NOERR) then
3658	write(,) NF_STRERROR(ierr)
3659	stop "getvarup"
3660	endif
3661	! write(,)'lecture qt ok',qt
3662
3663	#ifdef NC_DOUBLE
3664	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3665	#else
3666	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3667	#endif
3668	if(ierr/=NF_NOERR) then
3669	write(,) NF_STRERROR(ierr)
3670	stop "getvarup"
3671	endif
3672	! write(,)'lecture u ok',u
3673
3674	#ifdef NC_DOUBLE
3675	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3676	#else
3677	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3678	#endif
3679	if(ierr/=NF_NOERR) then
3680	write(,) NF_STRERROR(ierr)
3681	stop "getvarup"
3682	endif
3683	! write(,)'lecture v ok',v
3684
3685	#ifdef NC_DOUBLE
3686	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3687	#else
3688	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3689	#endif
3690	if(ierr/=NF_NOERR) then
3691	write(,) NF_STRERROR(ierr)
3692	stop "getvarup"
3693	endif
3694	! write(,)'lecture tke ok',tke
3695
3696	#ifdef NC_DOUBLE
3697	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3698	#else
3699	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3700	#endif
3701	if(ierr/=NF_NOERR) then
3702	write(,) NF_STRERROR(ierr)
3703	stop "getvarup"
3704	endif
3705	! write(,)'lecture ug ok',ug
3706
3707	#ifdef NC_DOUBLE
3708	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3709	#else
3710	ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3711	#endif
3712	if(ierr/=NF_NOERR) then
3713	write(,) NF_STRERROR(ierr)
3714	stop "getvarup"
3715	endif
3716	! write(,)'lecture vg ok',vg
3717
3718	#ifdef NC_DOUBLE
3719	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3720	#else
3721	ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3722	#endif
3723	if(ierr/=NF_NOERR) then
3724	write(,) NF_STRERROR(ierr)
3725	stop "getvarup"
3726	endif
3727	! write(,)'lecture wls ok',wls
3728
3729	#ifdef NC_DOUBLE
3730	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3731	#else
3732	ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3733	#endif
3734	if(ierr/=NF_NOERR) then
3735	write(,) NF_STRERROR(ierr)
3736	stop "getvarup"
3737	endif
3738	! write(,)'lecture dqtdx ok',dqtdx
3739
3740	#ifdef NC_DOUBLE
3741	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3742	#else
3743	ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3744	#endif
3745	if(ierr/=NF_NOERR) then
3746	write(,) NF_STRERROR(ierr)
3747	stop "getvarup"
3748	endif
3749	! write(,)'lecture dqtdy ok',dqtdy
3750
3751	#ifdef NC_DOUBLE
3752	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3753	#else
3754	ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3755	#endif
3756	if(ierr/=NF_NOERR) then
3757	write(,) NF_STRERROR(ierr)
3758	stop "getvarup"
3759	endif
3760	! write(,)'lecture dqtdt ok',dqtdt
3761
3762	#ifdef NC_DOUBLE
3763	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3764	#else
3765	ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3766	#endif
3767	if(ierr/=NF_NOERR) then
3768	write(,) NF_STRERROR(ierr)
3769	stop "getvarup"
3770	endif
3771	! write(,)'lecture thl_rad ok',thl_rad
3772
3773	return
3774	end subroutine read_fire

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

Download in other formats: