Context Navigation

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

Last change on this file since 1926 was 1914, checked in by Laurent Fairhead, 11 years ago

Modifs pour le cas FIRE
MPL

FIRE case mods
MPL

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.9 KB

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

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

Download in other formats: