Context Navigation

source: LMDZ5/branches/LF-private/libf/phy1d/1DUTILS.h @ 5405

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

La distinction avec ou sans writelim n'a plus lieu d'être
MP Lefebvre

the writelim/nowritelim distinction is of no further use
MP Lefebvre

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.

Download in other formats: