Context Navigation

source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/1DUTILS.h @ 2597

Last change on this file since 2597 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
File size: 127.6 KB

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

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

Download in other formats: