Context Navigation

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

Last change on this file since 1960 was 1960, checked in by fhourdin, 10 years ago

Nettoyage du code 1D pour limiter les warning en mode debug.
Essentiellement des declarations inutiles.

Cleaning of the 1D code in order to limit warnings in debug mode.
Essentially supressing unused declared variables.

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: 113.2 KB

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

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

Download in other formats: