Context Navigation

source: LMDZ5/branches/testing/libf/phylmd/1DUTILS.h @ 2212

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

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

Download in other formats: