Context Navigation

1DUTILS.h_no_writelim @ 2690

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

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL

AMMA and Sandu cases included in 1D configuration
MPL

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

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

Download in other formats:

Original Format