Context Navigation

1DUTILS.h @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 113.9 KB

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

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

Context Navigation

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

Download in other formats: