Context Navigation

source: LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 3241

Last change on this file since 3241 was 3008, checked in by Laurent Fairhead, 7 years ago

Modification in case theta exists in the forcings but not temperature.
We compute temp with theta

MPL

File size: 174.3 KB

Rev	Line
[2239]	1	#include "conf_gcm.F90"
[2018]	2
[2017]	3	!
	4	! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
	5	!
[2019]	6	!
	7	!
[2017]	8	SUBROUTINE conf_unicol
[2019]	9	!
[2017]	10	#ifdef CPP_IOIPSL
	11	use IOIPSL
	12	#else
	13	! if not using IOIPSL, we still need to use (a local version of) getin
	14	use ioipsl_getincom
	15	#endif
[2311]	16	USE print_control_mod, ONLY: lunout
[2017]	17	IMPLICIT NONE
[2019]	18	!-----------------------------------------------------------------------
	19	! Auteurs : A. Lahellec .
	20	!
	21	! Declarations :
	22	! --------------
[2017]	23
	24	#include "compar1d.h"
	25	#include "flux_arp.h"
	26	#include "tsoilnudge.h"
	27	#include "fcg_gcssold.h"
	28	#include "fcg_racmo.h"
[2019]	29	!
	30	!
	31	! local:
	32	! ------
[2017]	33
[2019]	34	! CHARACTER ch172,ch272,ch372,ch412
[2017]	35
[2019]	36	!
	37	! -------------------------------------------------------------------
	38	!
	39	! ......... Initilisation parametres du lmdz1D ..........
	40	!
	41	!---------------------------------------------------------------------
	42	! initialisations:
	43	! ----------------
[2017]	44
	45	!Config Key = lunout
	46	!Config Desc = unite de fichier pour les impressions
	47	!Config Def = 6
	48	!Config Help = unite de fichier pour les impressions
	49	!Config (defaut sortie standard = 6)
	50	lunout=6
	51	! CALL getin('lunout', lunout)
	52	IF (lunout /= 5 .and. lunout /= 6) THEN
	53	OPEN(lunout,FILE='lmdz.out')
	54	ENDIF
	55
	56	!Config Key = prt_level
[2716]	57	!Config Desc = niveau d'impressions de debogage
[2017]	58	!Config Def = 0
[2716]	59	!Config Help = Niveau d'impression pour le debogage
[2017]	60	!Config (0 = minimum d'impression)
[2019]	61	! prt_level = 0
	62	! CALL getin('prt_level',prt_level)
[2017]	63
[2019]	64	!-----------------------------------------------------------------------
	65	! Parametres de controle du run:
	66	!-----------------------------------------------------------------------
[2017]	67
	68	!Config Key = restart
	69	!Config Desc = on repart des startphy et start1dyn
	70	!Config Def = false
	71	!Config Help = les fichiers restart doivent etre renomme en start
[2019]	72	restart =.false.
[2017]	73	CALL getin('restart',restart)
	74
	75	!Config Key = forcing_type
	76	!Config Desc = defines the way the SCM is forced:
	77	!Config Def = 0
	78	!!Config Help = 0 ==> forcing_les = .true.
	79	! initial profiles from file prof.inp.001
	80	! no forcing by LS convergence ;
	81	! surface temperature imposed ;
	82	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
	83	! = 1 ==> forcing_radconv = .true.
	84	! idem forcing_type = 0, but the imposed radiative cooling
	85	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
	86	! then there is no radiative cooling at all)
	87	! = 2 ==> forcing_toga = .true.
	88	! initial profiles from TOGA-COARE IFA files
	89	! LS convergence and SST imposed from TOGA-COARE IFA files
	90	! = 3 ==> forcing_GCM2SCM = .true.
	91	! initial profiles from the GCM output
	92	! LS convergence imposed from the GCM output
	93	! = 4 ==> forcing_twpi = .true.
	94	! initial profiles from TWPICE nc files
	95	! LS convergence and SST imposed from TWPICE nc files
	96	! = 5 ==> forcing_rico = .true.
	97	! initial profiles from RICO idealized
	98	! LS convergence imposed from RICO (cst)
	99	! = 6 ==> forcing_amma = .true.
[2191]	100	! = 10 ==> forcing_case = .true.
	101	! initial profiles from case.nc file
[2017]	102	! = 40 ==> forcing_GCSSold = .true.
	103	! initial profile from GCSS file
	104	! LS convergence imposed from GCSS file
	105	! = 50 ==> forcing_fire = .true.
	106	! = 59 ==> forcing_sandu = .true.
	107	! initial profiles from sanduref file: see prof.inp.001
	108	! SST varying with time and divergence constante: see ifa_sanduref.txt file
	109	! Radiation has to be computed interactively
	110	! = 60 ==> forcing_astex = .true.
	111	! initial profiles from file: see prof.inp.001
	112	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
	113	! Radiation has to be computed interactively
	114	! = 61 ==> forcing_armcu = .true.
	115	! initial profiles from file: see prof.inp.001
	116	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
	117	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
	118	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
	119	! Radiation to be switched off
[2716]	120	! > 100 ==> forcing_case = .true. or forcing_case2 = .true.
	121	! initial profiles from case.nc file
[2017]	122	!
	123	forcing_type = 0
	124	CALL getin('forcing_type',forcing_type)
	125	imp_fcg_gcssold = .false.
	126	ts_fcg_gcssold = .false.
	127	Tp_fcg_gcssold = .false.
	128	Tp_ini_gcssold = .false.
	129	xTurb_fcg_gcssold = .false.
	130	IF (forcing_type .eq.40) THEN
	131	CALL getin('imp_fcg',imp_fcg_gcssold)
	132	CALL getin('ts_fcg',ts_fcg_gcssold)
	133	CALL getin('tp_fcg',Tp_fcg_gcssold)
	134	CALL getin('tp_ini',Tp_ini_gcssold)
	135	CALL getin('turb_fcg',xTurb_fcg_gcssold)
	136	ENDIF
	137
[2716]	138	!Parametres de forcage
[2191]	139	!Config Key = tend_t
	140	!Config Desc = forcage ou non par advection de T
	141	!Config Def = false
	142	!Config Help = forcage ou non par advection de T
	143	tend_t =0
	144	CALL getin('tend_t',tend_t)
	145
	146	!Config Key = tend_q
	147	!Config Desc = forcage ou non par advection de q
	148	!Config Def = false
	149	!Config Help = forcage ou non par advection de q
	150	tend_q =0
	151	CALL getin('tend_q',tend_q)
	152
	153	!Config Key = tend_u
	154	!Config Desc = forcage ou non par advection de u
	155	!Config Def = false
	156	!Config Help = forcage ou non par advection de u
	157	tend_u =0
	158	CALL getin('tend_u',tend_u)
	159
	160	!Config Key = tend_v
	161	!Config Desc = forcage ou non par advection de v
	162	!Config Def = false
	163	!Config Help = forcage ou non par advection de v
	164	tend_v =0
	165	CALL getin('tend_v',tend_v)
	166
	167	!Config Key = tend_w
	168	!Config Desc = forcage ou non par vitesse verticale
	169	!Config Def = false
	170	!Config Help = forcage ou non par vitesse verticale
	171	tend_w =0
	172	CALL getin('tend_w',tend_w)
	173
	174	!Config Key = tend_rayo
	175	!Config Desc = forcage ou non par dtrad
	176	!Config Def = false
	177	!Config Help = forcage ou non par dtrad
	178	tend_rayo =0
	179	CALL getin('tend_rayo',tend_rayo)
	180
	181
	182	!Config Key = nudge_t
	183	!Config Desc = constante de nudging de T
	184	!Config Def = false
	185	!Config Help = constante de nudging de T
	186	nudge_t =0.
	187	CALL getin('nudge_t',nudge_t)
	188
	189	!Config Key = nudge_q
	190	!Config Desc = constante de nudging de q
	191	!Config Def = false
	192	!Config Help = constante de nudging de q
	193	nudge_q =0.
	194	CALL getin('nudge_q',nudge_q)
	195
	196	!Config Key = nudge_u
	197	!Config Desc = constante de nudging de u
	198	!Config Def = false
	199	!Config Help = constante de nudging de u
	200	nudge_u =0.
	201	CALL getin('nudge_u',nudge_u)
	202
	203	!Config Key = nudge_v
	204	!Config Desc = constante de nudging de v
	205	!Config Def = false
	206	!Config Help = constante de nudging de v
	207	nudge_v =0.
	208	CALL getin('nudge_v',nudge_v)
	209
	210	!Config Key = nudge_w
	211	!Config Desc = constante de nudging de w
	212	!Config Def = false
	213	!Config Help = constante de nudging de w
	214	nudge_w =0.
	215	CALL getin('nudge_w',nudge_w)
	216
	217
[2181]	218	!Config Key = iflag_nudge
	219	!Config Desc = atmospheric nudging ttype (decimal code)
	220	!Config Def = 0
	221	!Config Help = 0 ==> no nudging
	222	! If digit number n of iflag_nudge is set, then nudging of type n is on
	223	! If digit number n of iflag_nudge is not set, then nudging of type n is off
	224	! (digits are numbered from the right)
	225	iflag_nudge = 0
	226	CALL getin('iflag_nudge',iflag_nudge)
	227
[2017]	228	!Config Key = ok_flux_surf
	229	!Config Desc = forcage ou non par les flux de surface
	230	!Config Def = false
	231	!Config Help = forcage ou non par les flux de surface
[2019]	232	ok_flux_surf =.false.
[2017]	233	CALL getin('ok_flux_surf',ok_flux_surf)
	234
[2126]	235	!Config Key = ok_prescr_ust
	236	!Config Desc = ustar impose ou non
	237	!Config Def = false
	238	!Config Help = ustar impose ou non
	239	ok_prescr_ust = .false.
	240	CALL getin('ok_prescr_ust',ok_prescr_ust)
	241
[2017]	242	!Config Key = ok_old_disvert
	243	!Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
	244	!Config Def = false
	245	!Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
	246	ok_old_disvert = .FALSE.
	247	CALL getin('ok_old_disvert',ok_old_disvert)
	248
	249	!Config Key = time_ini
	250	!Config Desc = meaningless in this case
	251	!Config Def = 0.
	252	!Config Help =
	253	tsurf = 0.
	254	CALL getin('time_ini',time_ini)
	255
	256	!Config Key = rlat et rlon
	257	!Config Desc = latitude et longitude
	258	!Config Def = 0.0 0.0
	259	!Config Help = fixe la position de la colonne
	260	xlat = 0.
	261	xlon = 0.
	262	CALL getin('rlat',xlat)
	263	CALL getin('rlon',xlon)
	264
	265	!Config Key = airephy
	266	!Config Desc = Grid cell area
	267	!Config Def = 1.e11
	268	!Config Help =
	269	airefi = 1.e11
	270	CALL getin('airephy',airefi)
	271
	272	!Config Key = nat_surf
	273	!Config Desc = surface type
	274	!Config Def = 0 (ocean)
	275	!Config Help = 0=ocean,1=land,2=glacier,3=banquise
	276	nat_surf = 0.
	277	CALL getin('nat_surf',nat_surf)
	278
	279	!Config Key = tsurf
	280	!Config Desc = surface temperature
	281	!Config Def = 290.
	282	!Config Help = not used if type_ts_forcing=1 in lmdz1d.F
	283	tsurf = 290.
	284	CALL getin('tsurf',tsurf)
	285
	286	!Config Key = psurf
	287	!Config Desc = surface pressure
	288	!Config Def = 102400.
	289	!Config Help =
	290	psurf = 102400.
	291	CALL getin('psurf',psurf)
	292
	293	!Config Key = zsurf
	294	!Config Desc = surface altitude
	295	!Config Def = 0.
	296	!Config Help =
	297	zsurf = 0.
	298	CALL getin('zsurf',zsurf)
	299
	300	!Config Key = rugos
	301	!Config Desc = coefficient de frottement
	302	!Config Def = 0.0001
	303	!Config Help = calcul du Cdrag
	304	rugos = 0.0001
	305	CALL getin('rugos',rugos)
	306
[2672]	307	!Config Key = rugosh
	308	!Config Desc = coefficient de frottement
	309	!Config Def = rugos
	310	!Config Help = calcul du Cdrag
	311	rugosh = rugos
	312	CALL getin('rugosh',rugosh)
	313
	314
	315
	316	!Config Key = snowmass
	317	!Config Desc = mass de neige de la surface en kg/m2
	318	!Config Def = 0.0000
	319	!Config Help = snowmass
	320	snowmass = 0.0000
	321	CALL getin('snowmass',snowmass)
	322
[2017]	323	!Config Key = wtsurf et wqsurf
	324	!Config Desc = ???
	325	!Config Def = 0.0 0.0
	326	!Config Help =
	327	wtsurf = 0.0
	328	wqsurf = 0.0
	329	CALL getin('wtsurf',wtsurf)
	330	CALL getin('wqsurf',wqsurf)
	331
	332	!Config Key = albedo
	333	!Config Desc = albedo
	334	!Config Def = 0.09
	335	!Config Help =
	336	albedo = 0.09
	337	CALL getin('albedo',albedo)
	338
	339	!Config Key = agesno
	340	!Config Desc = age de la neige
	341	!Config Def = 30.0
	342	!Config Help =
	343	xagesno = 30.0
	344	CALL getin('agesno',xagesno)
	345
	346	!Config Key = restart_runoff
	347	!Config Desc = age de la neige
	348	!Config Def = 30.0
	349	!Config Help =
	350	restart_runoff = 0.0
	351	CALL getin('restart_runoff',restart_runoff)
	352
	353	!Config Key = qsolinp
	354	!Config Desc = initial bucket water content (kg/m2) when land (5std)
	355	!Config Def = 30.0
	356	!Config Help =
	357	qsolinp = 1.
	358	CALL getin('qsolinp',qsolinp)
	359
	360	!Config Key = zpicinp
	361	!Config Desc = denivellation orographie
[2668]	362	!Config Def = 0.
[2017]	363	!Config Help = input brise
[2668]	364	zpicinp = 0.
[2017]	365	CALL getin('zpicinp',zpicinp)
	366	!Config key = nudge_tsoil
	367	!Config Desc = activation of soil temperature nudging
	368	!Config Def = .FALSE.
	369	!Config Help = ...
	370
	371	nudge_tsoil=.FALSE.
	372	CALL getin('nudge_tsoil',nudge_tsoil)
	373
	374	!Config key = isoil_nudge
	375	!Config Desc = level number where soil temperature is nudged
	376	!Config Def = 3
	377	!Config Help = ...
	378
	379	isoil_nudge=3
	380	CALL getin('isoil_nudge',isoil_nudge)
	381
	382	!Config key = Tsoil_nudge
	383	!Config Desc = target temperature for tsoil(isoil_nudge)
	384	!Config Def = 300.
	385	!Config Help = ...
	386
	387	Tsoil_nudge=300.
	388	CALL getin('Tsoil_nudge',Tsoil_nudge)
	389
	390	!Config key = tau_soil_nudge
	391	!Config Desc = nudging relaxation time for tsoil
	392	!Config Def = 3600.
	393	!Config Help = ...
	394
	395	tau_soil_nudge=3600.
	396	CALL getin('tau_soil_nudge',tau_soil_nudge)
	397
[2716]	398	!----------------------------------------------------------
	399	! Param??tres de for??age pour les forcages communs:
	400	! Pour les forcages communs: ces entiers valent 0 ou 1
	401	! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
	402	! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale
	403	! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv)
	404	! forcages en omega, w, vent geostrophique ou ustar
	405	! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging
	406	!----------------------------------------------------------
[2017]	407
[2716]	408	!Config Key = tadv
	409	!Config Desc = forcage ou non par advection totale de T
	410	!Config Def = false
	411	!Config Help = forcage ou non par advection totale de T
	412	tadv =0
	413	CALL getin('tadv',tadv)
[2017]	414
[2716]	415	!Config Key = tadvv
	416	!Config Desc = forcage ou non par advection verticale de T
	417	!Config Def = false
	418	!Config Help = forcage ou non par advection verticale de T
	419	tadvv =0
	420	CALL getin('tadvv',tadvv)
[2017]	421
[2716]	422	!Config Key = tadvh
	423	!Config Desc = forcage ou non par advection horizontale de T
	424	!Config Def = false
	425	!Config Help = forcage ou non par advection horizontale de T
	426	tadvh =0
	427	CALL getin('tadvh',tadvh)
	428
	429	!Config Key = thadv
	430	!Config Desc = forcage ou non par advection totale de Theta
	431	!Config Def = false
	432	!Config Help = forcage ou non par advection totale de Theta
	433	thadv =0
	434	CALL getin('thadv',thadv)
	435
	436	!Config Key = thadvv
	437	!Config Desc = forcage ou non par advection verticale de Theta
	438	!Config Def = false
	439	!Config Help = forcage ou non par advection verticale de Theta
	440	thadvv =0
	441	CALL getin('thadvv',thadvv)
	442
	443	!Config Key = thadvh
	444	!Config Desc = forcage ou non par advection horizontale de Theta
	445	!Config Def = false
	446	!Config Help = forcage ou non par advection horizontale de Theta
	447	thadvh =0
	448	CALL getin('thadvh',thadvh)
	449
	450	!Config Key = qadv
	451	!Config Desc = forcage ou non par advection totale de Q
	452	!Config Def = false
	453	!Config Help = forcage ou non par advection totale de Q
	454	qadv =0
	455	CALL getin('qadv',qadv)
	456
	457	!Config Key = qadvv
	458	!Config Desc = forcage ou non par advection verticale de Q
	459	!Config Def = false
	460	!Config Help = forcage ou non par advection verticale de Q
	461	qadvv =0
	462	CALL getin('qadvv',qadvv)
	463
	464	!Config Key = qadvh
	465	!Config Desc = forcage ou non par advection horizontale de Q
	466	!Config Def = false
	467	!Config Help = forcage ou non par advection horizontale de Q
	468	qadvh =0
	469	CALL getin('qadvh',qadvh)
	470
	471	!Config Key = trad
	472	!Config Desc = forcage ou non par tendance radiative
	473	!Config Def = false
	474	!Config Help = forcage ou non par tendance radiative
	475	trad =0
	476	CALL getin('trad',trad)
	477
	478	!Config Key = forc_omega
	479	!Config Desc = forcage ou non par omega
	480	!Config Def = false
	481	!Config Help = forcage ou non par omega
	482	forc_omega =0
	483	CALL getin('forc_omega',forc_omega)
	484
[2920]	485	!Config Key = forc_u
	486	!Config Desc = forcage ou non par u
	487	!Config Def = false
	488	!Config Help = forcage ou non par u
	489	forc_u =0
	490	CALL getin('forc_u',forc_u)
	491
	492	!Config Key = forc_v
	493	!Config Desc = forcage ou non par v
	494	!Config Def = false
	495	!Config Help = forcage ou non par v
	496	forc_v =0
	497	CALL getin('forc_v',forc_v)
[2716]	498	!Config Key = forc_w
	499	!Config Desc = forcage ou non par w
	500	!Config Def = false
	501	!Config Help = forcage ou non par w
	502	forc_w =0
	503	CALL getin('forc_w',forc_w)
	504
	505	!Config Key = forc_geo
	506	!Config Desc = forcage ou non par geo
	507	!Config Def = false
	508	!Config Help = forcage ou non par geo
	509	forc_geo =0
	510	CALL getin('forc_geo',forc_geo)
	511
	512	! Meme chose que ok_precr_ust
	513	!Config Key = forc_ustar
	514	!Config Desc = forcage ou non par ustar
	515	!Config Def = false
	516	!Config Help = forcage ou non par ustar
	517	forc_ustar =0
	518	CALL getin('forc_ustar',forc_ustar)
	519	IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
	520
	521	!Config Key = nudging_u
	522	!Config Desc = forcage ou non par nudging sur u
	523	!Config Def = false
	524	!Config Help = forcage ou non par nudging sur u
	525	nudging_u =0
	526	CALL getin('nudging_u',nudging_u)
	527
	528	!Config Key = nudging_v
	529	!Config Desc = forcage ou non par nudging sur v
	530	!Config Def = false
	531	!Config Help = forcage ou non par nudging sur v
	532	nudging_v =0
	533	CALL getin('nudging_v',nudging_v)
	534
	535	!Config Key = nudging_w
	536	!Config Desc = forcage ou non par nudging sur w
	537	!Config Def = false
	538	!Config Help = forcage ou non par nudging sur w
	539	nudging_w =0
	540	CALL getin('nudging_w',nudging_w)
	541
	542	!Config Key = nudging_q
	543	!Config Desc = forcage ou non par nudging sur q
	544	!Config Def = false
	545	!Config Help = forcage ou non par nudging sur q
	546	nudging_q =0
	547	CALL getin('nudging_q',nudging_q)
	548
	549	!Config Key = nudging_t
	550	!Config Desc = forcage ou non par nudging sur t
	551	!Config Def = false
	552	!Config Help = forcage ou non par nudging sur t
	553	nudging_t =0
	554	CALL getin('nudging_t',nudging_t)
	555
	556
	557
[2017]	558	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	559	write(lunout,*)' Configuration des parametres du gcm1D: '
	560	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	561	write(lunout,*)' restart = ', restart
	562	write(lunout,*)' forcing_type = ', forcing_type
	563	write(lunout,*)' time_ini = ', time_ini
	564	write(lunout,*)' rlat = ', xlat
	565	write(lunout,*)' rlon = ', xlon
	566	write(lunout,*)' airephy = ', airefi
	567	write(lunout,*)' nat_surf = ', nat_surf
	568	write(lunout,*)' tsurf = ', tsurf
	569	write(lunout,*)' psurf = ', psurf
	570	write(lunout,*)' zsurf = ', zsurf
	571	write(lunout,*)' rugos = ', rugos
[2672]	572	write(lunout,*)' snowmass=', snowmass
[2017]	573	write(lunout,*)' wtsurf = ', wtsurf
	574	write(lunout,*)' wqsurf = ', wqsurf
	575	write(lunout,*)' albedo = ', albedo
	576	write(lunout,*)' xagesno = ', xagesno
	577	write(lunout,*)' restart_runoff = ', restart_runoff
	578	write(lunout,*)' qsolinp = ', qsolinp
	579	write(lunout,*)' zpicinp = ', zpicinp
	580	write(lunout,*)' nudge_tsoil = ', nudge_tsoil
	581	write(lunout,*)' isoil_nudge = ', isoil_nudge
	582	write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
	583	write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
[2716]	584	write(lunout,*)' tadv = ', tadv
	585	write(lunout,*)' tadvv = ', tadvv
	586	write(lunout,*)' tadvh = ', tadvh
	587	write(lunout,*)' thadv = ', thadv
	588	write(lunout,*)' thadvv = ', thadvv
	589	write(lunout,*)' thadvh = ', thadvh
	590	write(lunout,*)' qadv = ', qadv
	591	write(lunout,*)' qadvv = ', qadvv
	592	write(lunout,*)' qadvh = ', qadvh
	593	write(lunout,*)' trad = ', trad
	594	write(lunout,*)' forc_omega = ', forc_omega
	595	write(lunout,*)' forc_w = ', forc_w
	596	write(lunout,*)' forc_geo = ', forc_geo
	597	write(lunout,*)' forc_ustar = ', forc_ustar
	598	write(lunout,*)' nudging_u = ', nudging_u
	599	write(lunout,*)' nudging_v = ', nudging_v
	600	write(lunout,*)' nudging_t = ', nudging_t
	601	write(lunout,*)' nudging_q = ', nudging_q
[2017]	602	IF (forcing_type .eq.40) THEN
	603	write(lunout,*) '--- Forcing type GCSS Old --- with:'
	604	write(lunout,*)'imp_fcg',imp_fcg_gcssold
	605	write(lunout,*)'ts_fcg',ts_fcg_gcssold
	606	write(lunout,*)'tp_fcg',Tp_fcg_gcssold
	607	write(lunout,*)'tp_ini',Tp_ini_gcssold
	608	write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
	609	ENDIF
	610
	611	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	612	write(lunout,*)
[2019]	613	!
[2017]	614	RETURN
	615	END
	616	!
	617	! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
	618	!
[2019]	619	!
	620	SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs, &
[2017]	621	& ucov,vcov,temp,q,omega2)
	622	USE dimphy
	623	USE mod_grid_phy_lmdz
	624	USE mod_phys_lmdz_para
	625	USE iophy
	626	USE phys_state_var_mod
	627	USE iostart
	628	USE write_field_phy
	629	USE infotrac
	630	use control_mod
[2597]	631	USE comconst_mod, ONLY: im, jm, lllm
[2603]	632	USE logic_mod, ONLY: fxyhypb, ysinus
[2601]	633	USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn
[2017]	634
	635	IMPLICIT NONE
[2019]	636	!=======================================================
	637	! Ecriture du fichier de redemarrage sous format NetCDF
	638	!=======================================================
	639	! Declarations:
	640	! -------------
[2239]	641	include "dimensions.h"
[2017]	642	!!#include "control.h"
[2239]	643	include "netcdf.inc"
[2017]	644
[2019]	645	! Arguments:
	646	! ----------
[2017]	647	CHARACTER() fichnom
[2019]	648	!Al1 plev tronque pour .nc mais plev(klev+1):=0
[2139]	649	real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev)
[2017]	650	real :: presnivs(klon,klev)
	651	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
	652	real :: q(klon,klev,nqtot),omega2(klon,klev)
[2019]	653	! real :: ug(klev),vg(klev),fcoriolis
[2017]	654	real :: phis(klon)
	655
[2019]	656	! Variables locales pour NetCDF:
	657	! ------------------------------
[2017]	658	INTEGER iq
	659	INTEGER length
	660	PARAMETER (length = 100)
	661	REAL tab_cntrl(length) ! tableau des parametres du run
	662	character*4 nmq(nqtot)
	663	character*12 modname
	664	character*80 abort_message
	665	LOGICAL found
	666
	667	modname = 'dyn1deta0 : '
[2933]	668	!! nmq(1)="vap"
	669	!! nmq(2)="cond"
	670	!! do iq=3,nqtot
	671	!! write(nmq(iq),'("tra",i1)') iq-2
	672	!! enddo
	673	DO iq = 1,nqtot
	674	nmq(iq) = trim(tname(iq))
	675	ENDDO
[2017]	676	print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
	677	CALL open_startphy(fichnom)
	678	print*,'after open startphy ',fichnom,nmq
	679
[2019]	680	!
	681	! Lecture des parametres de controle:
	682	!
[2017]	683	CALL get_var("controle",tab_cntrl)
	684
	685
	686	im = tab_cntrl(1)
	687	jm = tab_cntrl(2)
	688	lllm = tab_cntrl(3)
	689	day_ref = tab_cntrl(4)
	690	annee_ref = tab_cntrl(5)
[2019]	691	! rad = tab_cntrl(6)
	692	! omeg = tab_cntrl(7)
	693	! g = tab_cntrl(8)
	694	! cpp = tab_cntrl(9)
	695	! kappa = tab_cntrl(10)
	696	! daysec = tab_cntrl(11)
	697	! dtvr = tab_cntrl(12)
	698	! etot0 = tab_cntrl(13)
	699	! ptot0 = tab_cntrl(14)
	700	! ztot0 = tab_cntrl(15)
	701	! stot0 = tab_cntrl(16)
	702	! ang0 = tab_cntrl(17)
	703	! pa = tab_cntrl(18)
	704	! preff = tab_cntrl(19)
	705	!
	706	! clon = tab_cntrl(20)
	707	! clat = tab_cntrl(21)
	708	! grossismx = tab_cntrl(22)
	709	! grossismy = tab_cntrl(23)
	710	!
[2017]	711	IF ( tab_cntrl(24).EQ.1. ) THEN
[2019]	712	fxyhypb =.true.
	713	! dzoomx = tab_cntrl(25)
	714	! dzoomy = tab_cntrl(26)
	715	! taux = tab_cntrl(28)
	716	! tauy = tab_cntrl(29)
[2017]	717	ELSE
[2019]	718	fxyhypb = .false.
	719	ysinus = .false.
	720	IF( tab_cntrl(27).EQ.1. ) ysinus =.true.
[2017]	721	ENDIF
	722
	723	day_ini = tab_cntrl(30)
	724	itau_dyn = tab_cntrl(31)
[2019]	725	! .................................................................
	726	!
	727	!
	728	! PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
	729	!Al1
	730	Print*,'day_ref,annee_ref,day_ini,itau_dyn', &
[2017]	731	& day_ref,annee_ref,day_ini,itau_dyn
	732
[2019]	733	! Lecture des champs
	734	!
[2017]	735	CALL get_field("play",play,found)
	736	IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
	737	CALL get_field("phi",phi,found)
	738	IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
	739	CALL get_field("phis",phis,found)
	740	IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
	741	CALL get_field("presnivs",presnivs,found)
	742	IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
	743	CALL get_field("ucov",ucov,found)
	744	IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
	745	CALL get_field("vcov",vcov,found)
	746	IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
	747	CALL get_field("temp",temp,found)
	748	IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
	749	CALL get_field("omega2",omega2,found)
	750	IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
[2139]	751	plev(1,klev+1)=0.
	752	CALL get_field("plev",plev(:,1:klev),found)
	753	IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
[2017]	754
	755	Do iq=1,nqtot
	756	CALL get_field("q"//nmq(iq),q(:,:,iq),found)
[2019]	757	IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
[2017]	758	EndDo
	759
	760	CALL close_startphy
[2019]	761	print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
	762	!
[2017]	763	RETURN
	764	END
	765	!
	766	! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
	767	!
[2019]	768	!
	769	SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs, &
[2017]	770	& ucov,vcov,temp,q,omega2)
	771	USE dimphy
	772	USE mod_grid_phy_lmdz
	773	USE mod_phys_lmdz_para
	774	USE phys_state_var_mod
	775	USE iostart
	776	USE infotrac
	777	use control_mod
[2597]	778	USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
[2603]	779	USE logic_mod, ONLY: fxyhypb, ysinus
[2601]	780	USE temps_mod, ONLY: annee_ref,day_end,day_ref,itau_dyn,itaufin
[2017]	781
	782	IMPLICIT NONE
[2019]	783	!=======================================================
	784	! Ecriture du fichier de redemarrage sous format NetCDF
	785	!=======================================================
	786	! Declarations:
	787	! -------------
[2239]	788	include "dimensions.h"
[2017]	789	!!#include "control.h"
[2239]	790	include "netcdf.inc"
[2017]	791
[2019]	792	! Arguments:
	793	! ----------
[2017]	794	CHARACTER() fichnom
[2019]	795	!Al1 plev tronque pour .nc mais plev(klev+1):=0
[2017]	796	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
	797	real :: presnivs(klon,klev)
	798	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
	799	real :: q(klon,klev,nqtot)
	800	real :: omega2(klon,klev),rho(klon,klev+1)
[2019]	801	! real :: ug(klev),vg(klev),fcoriolis
[2017]	802	real :: phis(klon)
	803
[2019]	804	! Variables locales pour NetCDF:
	805	! ------------------------------
[2017]	806	INTEGER nid
	807	INTEGER ierr
	808	INTEGER iq,l
	809	INTEGER length
	810	PARAMETER (length = 100)
	811	REAL tab_cntrl(length) ! tableau des parametres du run
	812	character*4 nmq(nqtot)
	813	character*20 modname
	814	character*80 abort_message
[2019]	815	!
[2017]	816	INTEGER nb
	817	SAVE nb
	818	DATA nb / 0 /
	819
	820	CALL open_restartphy(fichnom)
	821	print*,'redm1 ',fichnom,klon,klev,nqtot
[2933]	822	!! nmq(1)="vap"
	823	!! nmq(2)="cond"
	824	!! nmq(3)="tra1"
	825	!! nmq(4)="tra2"
	826	DO iq = 1,nqtot
	827	nmq(iq) = trim(tname(iq))
	828	ENDDO
[2017]	829
	830	modname = 'dyn1dredem'
	831	ierr = NF_OPEN(fichnom, NF_WRITE, nid)
	832	IF (ierr .NE. NF_NOERR) THEN
	833	abort_message="Pb. d ouverture "//fichnom
	834	CALL abort_gcm('Modele 1D',abort_message,1)
	835	ENDIF
	836
	837	DO l=1,length
	838	tab_cntrl(l) = 0.
	839	ENDDO
	840	tab_cntrl(1) = FLOAT(iim)
	841	tab_cntrl(2) = FLOAT(jjm)
	842	tab_cntrl(3) = FLOAT(llm)
	843	tab_cntrl(4) = FLOAT(day_ref)
	844	tab_cntrl(5) = FLOAT(annee_ref)
	845	tab_cntrl(6) = rad
	846	tab_cntrl(7) = omeg
	847	tab_cntrl(8) = g
	848	tab_cntrl(9) = cpp
	849	tab_cntrl(10) = kappa
	850	tab_cntrl(11) = daysec
	851	tab_cntrl(12) = dtvr
[2019]	852	! tab_cntrl(13) = etot0
	853	! tab_cntrl(14) = ptot0
	854	! tab_cntrl(15) = ztot0
	855	! tab_cntrl(16) = stot0
	856	! tab_cntrl(17) = ang0
	857	! tab_cntrl(18) = pa
	858	! tab_cntrl(19) = preff
	859	!
	860	! ..... parametres pour le zoom ......
[2017]	861
[2019]	862	! tab_cntrl(20) = clon
	863	! tab_cntrl(21) = clat
	864	! tab_cntrl(22) = grossismx
	865	! tab_cntrl(23) = grossismy
	866	!
[2017]	867	IF ( fxyhypb ) THEN
	868	tab_cntrl(24) = 1.
[2019]	869	! tab_cntrl(25) = dzoomx
	870	! tab_cntrl(26) = dzoomy
[2017]	871	tab_cntrl(27) = 0.
[2019]	872	! tab_cntrl(28) = taux
	873	! tab_cntrl(29) = tauy
[2017]	874	ELSE
	875	tab_cntrl(24) = 0.
[2019]	876	! tab_cntrl(25) = dzoomx
	877	! tab_cntrl(26) = dzoomy
[2017]	878	tab_cntrl(27) = 0.
	879	tab_cntrl(28) = 0.
	880	tab_cntrl(29) = 0.
	881	IF( ysinus ) tab_cntrl(27) = 1.
	882	ENDIF
[2019]	883	!Al1 iday_end -> day_end
[2017]	884	tab_cntrl(30) = FLOAT(day_end)
	885	tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
[2019]	886	!
[2017]	887	CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
[2019]	888	!
[2017]	889
[2019]	890	! Ecriture/extension de la coordonnee temps
[2017]	891
	892	nb = nb + 1
	893
[2019]	894	! Ecriture des champs
	895	!
[2017]	896	CALL put_field("plev","p interfaces sauf la nulle",plev)
	897	CALL put_field("play","",play)
	898	CALL put_field("phi","geopotentielle",phi)
	899	CALL put_field("phis","geopotentiell de surface",phis)
	900	CALL put_field("presnivs","",presnivs)
	901	CALL put_field("ucov","",ucov)
	902	CALL put_field("vcov","",vcov)
	903	CALL put_field("temp","",temp)
	904	CALL put_field("omega2","",omega2)
	905
	906	Do iq=1,nqtot
[2019]	907	CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs", &
	908	& q(:,:,iq))
[2017]	909	EndDo
	910	CALL close_restartphy
	911
[2019]	912	!
[2017]	913	RETURN
	914	END
	915	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
	916	IMPLICIT NONE
	917	!=======================================================================
	918	! passage d'un champ de la grille scalaire a la grille physique
	919	!=======================================================================
	920
	921	!-----------------------------------------------------------------------
	922	! declarations:
	923	! -------------
	924
	925	INTEGER im,jm,ngrid,nfield
	926	REAL pdyn(im,jm,nfield)
	927	REAL pfi(ngrid,nfield)
	928
	929	INTEGER i,j,ifield,ig
	930
	931	!-----------------------------------------------------------------------
	932	! calcul:
	933	! -------
	934
	935	DO ifield=1,nfield
	936	! traitement des poles
	937	DO i=1,im
	938	pdyn(i,1,ifield)=pfi(1,ifield)
	939	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
	940	ENDDO
	941
	942	! traitement des point normaux
	943	DO j=2,jm-1
	944	ig=2+(j-2)*(im-1)
	945	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
	946	pdyn(im,j,ifield)=pdyn(1,j,ifield)
	947	ENDDO
	948	ENDDO
	949
	950	RETURN
	951	END
	952
	953
	954
	955	SUBROUTINE abort_gcm(modname, message, ierr)
	956
	957	USE IOIPSL
	958	!
	959	! Stops the simulation cleanly, closing files and printing various
	960	! comments
	961	!
	962	! Input: modname = name of calling program
	963	! message = stuff to print
	964	! ierr = severity of situation ( = 0 normal )
	965
	966	character(len=*) modname
	967	integer ierr
	968	character(len=*) message
	969
	970	write(,) 'in abort_gcm'
	971	call histclo
	972	! call histclo(2)
	973	! call histclo(3)
	974	! call histclo(4)
	975	! call histclo(5)
	976	write(,) 'out of histclo'
	977	write(,) 'Stopping in ', modname
	978	write(,) 'Reason = ',message
	979	call getin_dump
	980	!
	981	if (ierr .eq. 0) then
	982	write(,) 'Everything is cool'
	983	else
	984	write(,) 'Houston, we have a problem ', ierr
	985	endif
	986	STOP
	987	END
	988	REAL FUNCTION fq_sat(kelvin, millibar)
	989	!
	990	IMPLICIT none
	991	!======================================================================
	992	! Autheur(s): Z.X. Li (LMD/CNRS)
	993	! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
	994	!======================================================================
	995	! Arguments:
	996	! kelvin---input-R: temperature en Kelvin
	997	! millibar--input-R: pression en mb
	998	!
	999	! fq_sat----output-R: vapeur d'eau saturante en kg/kg
	1000	!======================================================================
	1001	!
	1002	REAL kelvin, millibar
	1003	!
	1004	REAL r2es
	1005	PARAMETER (r2es=611.14 *18.0153/28.9644)
	1006	!
	1007	REAL r3les, r3ies, r3es
	1008	PARAMETER (R3LES=17.269)
	1009	PARAMETER (R3IES=21.875)
	1010	!
	1011	REAL r4les, r4ies, r4es
	1012	PARAMETER (R4LES=35.86)
	1013	PARAMETER (R4IES=7.66)
	1014	!
	1015	REAL rtt
	1016	PARAMETER (rtt=273.16)
	1017	!
	1018	REAL retv
	1019	PARAMETER (retv=28.9644/18.0153 - 1.0)
	1020	!
	1021	REAL zqsat
	1022	REAL temp, pres
	1023	! ------------------------------------------------------------------
	1024	!
	1025	!
	1026	temp = kelvin
	1027	pres = millibar * 100.0
	1028	! write(,)'kelvin,millibar=',kelvin,millibar
	1029	! write(,)'temp,pres=',temp,pres
	1030	!
	1031	IF (temp .LE. rtt) THEN
	1032	r3es = r3ies
	1033	r4es = r4ies
	1034	ELSE
	1035	r3es = r3les
	1036	r4es = r4les
	1037	ENDIF
	1038	!
	1039	zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
	1040	zqsat=MIN(0.5,ZQSAT)
	1041	zqsat=zqsat/(1.-retv *zqsat)
	1042	!
	1043	fq_sat = zqsat
	1044	!
	1045	RETURN
	1046	END
	1047
	1048	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
	1049	IMPLICIT NONE
	1050	!=======================================================================
	1051	! passage d'un champ de la grille scalaire a la grille physique
	1052	!=======================================================================
	1053
	1054	!-----------------------------------------------------------------------
	1055	! declarations:
	1056	! -------------
	1057
	1058	INTEGER im,jm,ngrid,nfield
	1059	REAL pdyn(im,jm,nfield)
	1060	REAL pfi(ngrid,nfield)
	1061
	1062	INTEGER j,ifield,ig
	1063
	1064	!-----------------------------------------------------------------------
	1065	! calcul:
	1066	! -------
	1067
[2019]	1068	IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1) &
	1069	& STOP 'probleme de dim'
[2017]	1070	! traitement des poles
	1071	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
	1072	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
	1073
	1074	! traitement des point normaux
	1075	DO ifield=1,nfield
	1076	DO j=2,jm-1
	1077	ig=2+(j-2)*(im-1)
	1078	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
	1079	ENDDO
	1080	ENDDO
	1081
	1082	RETURN
	1083	END
	1084
	1085	SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
	1086
	1087	! Ancienne version disvert dont on a modifie nom pour utiliser
	1088	! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
	1089	! (MPL 18092012)
	1090	!
	1091	! Auteur : P. Le Van .
	1092	!
	1093	IMPLICIT NONE
	1094
[2239]	1095	include "dimensions.h"
	1096	include "paramet.h"
[2017]	1097	!
	1098	!=======================================================================
	1099	!
	1100	!
	1101	! s = sigma ** kappa : coordonnee verticale
	1102	! dsig(l) : epaisseur de la couche l ds la coord. s
	1103	! sig(l) : sigma a l'interface des couches l et l-1
	1104	! ds(l) : distance entre les couches l et l-1 en coord.s
	1105	!
	1106	!=======================================================================
	1107	!
	1108	REAL pa,preff
	1109	REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
	1110	REAL presnivs(llm)
	1111	!
	1112	! declarations:
	1113	! -------------
	1114	!
	1115	REAL sig(llm+1),dsig(llm)
	1116	!
	1117	INTEGER l
	1118	REAL snorm
	1119	REAL alpha,beta,gama,delta,deltaz,h
	1120	INTEGER np,ierr
	1121	REAL pi,x
	1122
	1123	!-----------------------------------------------------------------------
	1124	!
	1125	pi=2.*ASIN(1.)
	1126
[2019]	1127	OPEN(99,file='sigma.def',status='old',form='formatted', &
	1128	& iostat=ierr)
[2017]	1129
	1130	!-----------------------------------------------------------------------
	1131	! cas 1 on lit les options dans sigma.def:
	1132	! ----------------------------------------
	1133
	1134	IF (ierr.eq.0) THEN
	1135
	1136	print*,'WARNING!!! on lit les options dans sigma.def'
	1137	READ(99,*) deltaz
	1138	READ(99,*) h
	1139	READ(99,*) beta
	1140	READ(99,*) gama
	1141	READ(99,*) delta
	1142	READ(99,*) np
	1143	CLOSE(99)
	1144	alpha=deltaz/(llm*h)
	1145	!
	1146
	1147	DO 1 l = 1, llm
[2019]	1148	dsig(l) = (alpha+(1.-alpha)exp(-beta(llm-l)))* &
	1149	& ( (tanh(gamal)/tanh(gamallm))**np + &
	1150	& (1.-l/FLOAT(llm))*delta )
[2017]	1151	1 CONTINUE
	1152
	1153	sig(1)=1.
	1154	DO 101 l=1,llm-1
	1155	sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
	1156	101 CONTINUE
	1157	sig(llm+1)=0.
	1158
	1159	DO 2 l = 1, llm
	1160	dsig(l) = sig(l)-sig(l+1)
	1161	2 CONTINUE
	1162	!
	1163
	1164	ELSE
	1165	!-----------------------------------------------------------------------
	1166	! cas 2 ancienne discretisation (LMD5...):
	1167	! ----------------------------------------
	1168
	1169	PRINT*,'WARNING!!! Ancienne discretisation verticale'
	1170
	1171	h=7.
	1172	snorm = 0.
	1173	DO l = 1, llm
	1174	x = 2.asin(1.) (FLOAT(l)-0.5) / float(llm+1)
	1175	dsig(l) = 1.0 + 7.0 * SIN(x)**2
	1176	snorm = snorm + dsig(l)
	1177	ENDDO
	1178	snorm = 1./snorm
	1179	DO l = 1, llm
	1180	dsig(l) = dsig(l)*snorm
	1181	ENDDO
	1182	sig(llm+1) = 0.
	1183	DO l = llm, 1, -1
	1184	sig(l) = sig(l+1) + dsig(l)
	1185	ENDDO
	1186
	1187	ENDIF
	1188
	1189
	1190	DO l=1,llm
	1191	nivsigs(l) = FLOAT(l)
	1192	ENDDO
	1193
	1194	DO l=1,llmp1
	1195	nivsig(l)= FLOAT(l)
	1196	ENDDO
	1197
	1198	!
	1199	! .... Calculs de ap(l) et de bp(l) ....
	1200	! .........................................
	1201	!
	1202	!
	1203	! ..... pa et preff sont lus sur les fichiers start par lectba .....
	1204	!
	1205
	1206	bp(llmp1) = 0.
	1207
	1208	DO l = 1, llm
	1209	!c
	1210	!cc ap(l) = 0.
	1211	!cc bp(l) = sig(l)
	1212
	1213	bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
	1214	ap(l) = pa * ( sig(l) - bp(l) )
	1215	!
	1216	ENDDO
	1217	ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
	1218
	1219	PRINT *,' BP '
	1220	PRINT *, bp
	1221	PRINT *,' AP '
	1222	PRINT *, ap
	1223
[2096]	1224	DO l = 1, llm
	1225	dpres(l) = bp(l) - bp(l+1)
	1226	presnivs(l) = 0.5 ( ap(l)+bp(l)preff + ap(l+1)+bp(l+1)*preff )
	1227	ENDDO
[2017]	1228
	1229	PRINT *,' PRESNIVS '
	1230	PRINT *,presnivs
	1231
	1232	RETURN
	1233	END
	1234
	1235	!======================================================================
	1236	SUBROUTINE read_tsurf1d(knon,sst_out)
	1237
	1238	! This subroutine specifies the surface temperature to be used in 1D simulations
	1239
	1240	USE dimphy, ONLY : klon
	1241
	1242	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
	1243	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
	1244
	1245	INTEGER :: i
	1246	! COMMON defined in lmdz1d.F:
	1247	real ts_cur
	1248	common /sst_forcing/ts_cur
	1249
	1250	DO i = 1, knon
	1251	sst_out(i) = ts_cur
	1252	ENDDO
	1253
	1254	END SUBROUTINE read_tsurf1d
	1255
	1256	!===============================================================
	1257	subroutine advect_vert(llm,w,dt,q,plev)
	1258	!===============================================================
	1259	! Schema amont pour l'advection verticale en 1D
	1260	! w est la vitesse verticale dp/dt en Pa/s
	1261	! Traitement en volumes finis
	1262	! d / dt ( zm q ) = delta_z ( omega q )
	1263	! d / dt ( zm ) = delta_z ( omega )
	1264	! avec zm = delta_z ( p )
	1265	! si * designe la valeur au pas de temps t+dt
	1266	! zm(l) q(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
	1267	! zm*(l) -zm(l) = w(l+1) - w(l)
	1268	! avec w=omega * dt
	1269	!---------------------------------------------------------------
	1270	implicit none
	1271	! arguments
	1272	integer llm
	1273	real w(llm+1),q(llm),plev(llm+1),dt
	1274
	1275	! local
	1276	integer l
	1277	real zwq(llm+1),zm(llm+1),zw(llm+1)
	1278	real qold
	1279
	1280	!---------------------------------------------------------------
	1281
	1282	do l=1,llm
	1283	zw(l)=dt*w(l)
	1284	zm(l)=plev(l)-plev(l+1)
	1285	zwq(l)=q(l)*zw(l)
	1286	enddo
	1287	zwq(llm+1)=0.
	1288	zw(llm+1)=0.
	1289
	1290	do l=1,llm
	1291	qold=q(l)
	1292	q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
	1293	print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
	1294	enddo
	1295
	1296
	1297	return
	1298	end
	1299
	1300	!===============================================================
	1301
	1302
[2019]	1303	SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, &
	1304	& q,temp,u,v,play)
[2017]	1305	!itlmd
	1306	!----------------------------------------------------------------------
	1307	! Calcul de l'advection verticale (ascendance et subsidence) de
[2716]	1308	! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
	1309	! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
[2017]	1310	! sans WTG rajouter une advection horizontale
	1311	!----------------------------------------------------------------------
	1312	implicit none
	1313	#include "YOMCST.h"
	1314	! argument
	1315	integer llm
	1316	real omega(llm+1),d_t_va(llm), d_q_va(llm,3)
	1317	real d_u_va(llm), d_v_va(llm)
	1318	real q(llm,3),temp(llm)
	1319	real u(llm),v(llm)
	1320	real play(llm)
	1321	! interne
	1322	integer l
	1323	real alpha,omgdown,omgup
	1324
	1325	do l= 1,llm
	1326	if(l.eq.1) then
	1327	!si omgup pour la couche 1, alors tendance nulle
	1328	omgdown=max(omega(2),0.0)
[2019]	1329	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
	1330	d_t_va(l)= alpha(omgdown)-omgdown(temp(l)-temp(l+1)) &
	1331	& /(play(l)-play(l+1))
[2017]	1332
[2019]	1333	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))
[2017]	1334
[2019]	1335	d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))
	1336	d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))
[2017]	1337
	1338
	1339	elseif(l.eq.llm) then
	1340	omgup=min(omega(l),0.0)
[2019]	1341	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
	1342	d_t_va(l)= alpha*(omgup)- &
	1343
[2017]	1344	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
	1345	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
	1346	d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
	1347	d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
	1348	d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
	1349
	1350	else
	1351	omgup=min(omega(l),0.0)
	1352	omgdown=max(omega(l+1),0.0)
[2019]	1353	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)(1.+q(l,1)))
	1354	d_t_va(l)= alpha(omgup+omgdown)-omgdown(temp(l)-temp(l+1)) &
	1355	& /(play(l)-play(l+1))- &
[2017]	1356	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
	1357	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
	1358	! print*, ' ??? '
	1359
[2019]	1360	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:)) &
	1361	& /(play(l)-play(l+1))- &
	1362	& omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
	1363	d_u_va(l)= -omgdown*(u(l)-u(l+1)) &
	1364	& /(play(l)-play(l+1))- &
	1365	& omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
	1366	d_v_va(l)= -omgdown*(v(l)-v(l+1)) &
	1367	& /(play(l)-play(l+1))- &
[2017]	1368	& omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
	1369
	1370	endif
	1371
	1372	enddo
	1373	!fin itlmd
	1374	return
	1375	end
	1376	! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
[2019]	1377	SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va, &
	1378	& q,temp,u,v,play)
[2017]	1379	!itlmd
	1380	!----------------------------------------------------------------------
	1381	! Calcul de l'advection verticale (ascendance et subsidence) de
[2716]	1382	! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
	1383	! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
[2017]	1384	! sans WTG rajouter une advection horizontale
	1385	!----------------------------------------------------------------------
	1386	implicit none
	1387	#include "YOMCST.h"
	1388	! argument
	1389	integer llm,nqtot
	1390	real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
	1391	! real d_u_va(llm), d_v_va(llm)
	1392	real q(llm,nqtot),temp(llm)
	1393	real u(llm),v(llm)
	1394	real play(llm)
	1395	real cor(llm)
	1396	! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
	1397	real dph(llm),dqdp(llm),dtdp(llm)
	1398	! interne
	1399	integer k
	1400	real omdn,omup
	1401
	1402	! dudp=0.
	1403	! dvdp=0.
	1404	dqdp=0.
	1405	dtdp=0.
	1406	! d_u_va=0.
	1407	! d_v_va=0.
	1408
	1409	cor(:) = rkappatemp(1.+q(:,1)rv/rd)/(play(1.+q(:,1)))
	1410
[2933]	1411
[2017]	1412	do k=2,llm-1
	1413
	1414	dph (k-1) = (play(k )- play(k-1 ))
	1415	! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
	1416	! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1)
	1417	dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1)
	1418	dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1)
	1419
	1420	enddo
	1421
	1422	! dudp ( llm ) = dudp ( llm-1 )
	1423	! dvdp ( llm ) = dvdp ( llm-1 )
	1424	dqdp ( llm ) = dqdp ( llm-1 )
	1425	dtdp ( llm ) = dtdp ( llm-1 )
	1426
	1427	do k=2,llm-1
	1428	omdn=max(0.0,omega(k+1))
	1429	omup=min(0.0,omega( k ))
	1430
	1431	! d_u_va(k) = -omdndudp(k)-omupdudp(k-1)
	1432	! d_v_va(k) = -omdndvdp(k)-omupdvdp(k-1)
	1433	d_q_va(k,1)= -omdndqdp(k)-omupdqdp(k-1)
[2019]	1434	d_t_va(k) = -omdndtdp(k)-omupdtdp(k-1)+(omup+omdn)*cor(k)
[2017]	1435	enddo
	1436
	1437	omdn=max(0.0,omega( 2 ))
	1438	omup=min(0.0,omega(llm))
	1439	! d_u_va( 1 ) = -omdn*dudp( 1 )
	1440	! d_u_va(llm) = -omup*dudp(llm)
	1441	! d_v_va( 1 ) = -omdn*dvdp( 1 )
	1442	! d_v_va(llm) = -omup*dvdp(llm)
	1443	d_q_va( 1 ,1) = -omdn*dqdp( 1 )
	1444	d_q_va(llm,1) = -omup*dqdp(llm)
	1445	d_t_va( 1 ) = -omdndtdp( 1 )+omdncor( 1 )
	1446	d_t_va(llm) = -omupdtdp(llm)!+omupcor(llm)
	1447
	1448	! if(abs(rlat(1))>10.) then
	1449	! Calculate the tendency due agestrophic motions
	1450	! du_age = fcoriolis*(v-vg)
	1451	! dv_age = fcoriolis*(ug-u)
	1452	! endif
	1453
	1454	! call writefield_phy('d_t_va',d_t_va,llm)
	1455
	1456	return
	1457	end
	1458
	1459	!======================================================================
[2019]	1460	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &
	1461	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
	1462	& ,ht_toga,vt_toga,hq_toga,vq_toga)
[2017]	1463	implicit none
	1464
[2019]	1465	!-------------------------------------------------------------------------
	1466	! Read TOGA-COARE forcing data
	1467	!-------------------------------------------------------------------------
[2017]	1468
	1469	integer nlev_toga,nt_toga
	1470	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
	1471	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
	1472	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
	1473	real w_toga(nlev_toga,nt_toga)
	1474	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
	1475	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
	1476	character*80 fich_toga
	1477
	1478	integer k,ip
	1479	real bid
	1480
	1481	integer iy,im,id,ih
	1482
	1483	real plev_min
	1484
	1485	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
	1486
	1487	open(21,file=trim(fich_toga),form='formatted')
	1488	read(21,'(a)')
	1489	do ip = 1, nt_toga
	1490	read(21,'(a)')
	1491	read(21,'(a)')
	1492	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
	1493	read(21,'(a)')
	1494	read(21,'(a)')
	1495
	1496	do k = 1, nlev_toga
[2019]	1497	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &
	1498	& ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &
	1499	& ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
[2017]	1500
	1501	! conversion in SI units:
	1502	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
	1503	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
	1504	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
	1505	! no water vapour tendency above 55 hPa
	1506	if (plev_toga(k,ip) .lt. plev_min) then
	1507	q_toga(k,ip) = 0.
	1508	hq_toga(k,ip) = 0.
	1509	vq_toga(k,ip) =0.
	1510	endif
	1511	enddo
	1512
	1513	ts_toga(ip)=ts_toga(ip)+273.15 ! K
	1514	enddo
	1515	close(21)
	1516
	1517	223 format(4i3,6f8.2)
	1518	230 format(6f9.3,4e11.3)
	1519
	1520	return
	1521	end
	1522
[2019]	1523	!-------------------------------------------------------------------------
[2017]	1524	SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
	1525	implicit none
	1526
[2019]	1527	!-------------------------------------------------------------------------
	1528	! Read I.SANDU case forcing data
	1529	!-------------------------------------------------------------------------
[2017]	1530
	1531	integer nlev_sandu,nt_sandu
	1532	real ts_sandu(nt_sandu)
	1533	character*80 fich_sandu
	1534
	1535	integer ip
	1536	integer iy,im,id,ih
	1537
	1538	real plev_min
	1539
	1540	print*,'nlev_sandu',nlev_sandu
	1541	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
	1542
	1543	open(21,file=trim(fich_sandu),form='formatted')
	1544	read(21,'(a)')
	1545	do ip = 1, nt_sandu
	1546	read(21,'(a)')
	1547	read(21,'(a)')
	1548	read(21,223) iy, im, id, ih, ts_sandu(ip)
	1549	print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
	1550	enddo
	1551	close(21)
	1552
	1553	223 format(4i3,f8.2)
	1554
	1555	return
	1556	end
	1557
	1558	!=====================================================================
[2019]	1559	!-------------------------------------------------------------------------
	1560	SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, &
	1561	& ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
[2017]	1562	implicit none
	1563
[2019]	1564	!-------------------------------------------------------------------------
	1565	! Read Astex case forcing data
	1566	!-------------------------------------------------------------------------
[2017]	1567
	1568	integer nlev_astex,nt_astex
	1569	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
	1570	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
	1571	character*80 fich_astex
	1572
	1573	integer ip
	1574	integer iy,im,id,ih
	1575
	1576	real plev_min
	1577
	1578	print*,'nlev_astex',nlev_astex
	1579	plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa
	1580
	1581	open(21,file=trim(fich_astex),form='formatted')
	1582	read(21,'(a)')
	1583	read(21,'(a)')
	1584	do ip = 1, nt_astex
	1585	read(21,'(a)')
	1586	read(21,'(a)')
[2019]	1587	read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), &
	1588	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
[2017]	1589	ts_astex(ip)=ts_astex(ip)+273.15
[2019]	1590	print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), &
	1591	&ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
[2017]	1592	enddo
	1593	close(21)
	1594
	1595	223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
	1596
	1597	return
	1598	end
	1599	!=====================================================================
[2019]	1600	subroutine read_twpice(fich_twpice,nlevel,ntime &
	1601	& ,T_srf,plev,T,q,u,v,omega &
	1602	& ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
[2017]	1603
	1604	!program reading forcings of the TWP-ICE experiment
	1605
	1606	! use netcdf
	1607
	1608	implicit none
	1609
	1610	#include "netcdf.inc"
	1611
	1612	integer ntime,nlevel
	1613	integer l,k
	1614	character*80 :: fich_twpice
	1615	real*8 time(ntime)
	1616	real*8 lat, lon, alt, phis
	1617	real*8 lev(nlevel)
	1618	real*8 plev(nlevel,ntime)
	1619
	1620	real*8 T(nlevel,ntime)
	1621	real*8 q(nlevel,ntime),u(nlevel,ntime)
	1622	real*8 v(nlevel,ntime)
	1623	real*8 omega(nlevel,ntime), div(nlevel,ntime)
	1624	real*8 T_adv_h(nlevel,ntime)
	1625	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
	1626	real*8 q_adv_v(nlevel,ntime)
	1627	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
	1628	real*8 s_adv_v(nlevel,ntime)
	1629	real*8 p_srf_aver(ntime), p_srf_center(ntime)
	1630	real*8 T_srf(ntime)
	1631
	1632	integer nid, ierr
	1633	integer nbvar3d
	1634	parameter(nbvar3d=20)
	1635	integer var3didin(nbvar3d)
	1636
	1637	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
	1638	if (ierr.NE.NF_NOERR) then
	1639	write(,) 'ERROR: Pb opening forcings cdf file '
	1640	write(,) NF_STRERROR(ierr)
	1641	stop ""
	1642	endif
	1643
	1644	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
	1645	if(ierr/=NF_NOERR) then
	1646	write(,) NF_STRERROR(ierr)
	1647	stop 'lat'
	1648	endif
	1649
	1650	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
	1651	if(ierr/=NF_NOERR) then
	1652	write(,) NF_STRERROR(ierr)
	1653	stop 'lon'
	1654	endif
	1655
	1656	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
	1657	if(ierr/=NF_NOERR) then
	1658	write(,) NF_STRERROR(ierr)
	1659	stop 'alt'
	1660	endif
	1661
	1662	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
	1663	if(ierr/=NF_NOERR) then
	1664	write(,) NF_STRERROR(ierr)
	1665	stop 'phis'
	1666	endif
	1667
	1668	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
	1669	if(ierr/=NF_NOERR) then
	1670	write(,) NF_STRERROR(ierr)
	1671	stop 'T'
	1672	endif
	1673
	1674	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
	1675	if(ierr/=NF_NOERR) then
	1676	write(,) NF_STRERROR(ierr)
	1677	stop 'q'
	1678	endif
	1679
	1680	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
	1681	if(ierr/=NF_NOERR) then
	1682	write(,) NF_STRERROR(ierr)
	1683	stop 'u'
	1684	endif
	1685
	1686	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
	1687	if(ierr/=NF_NOERR) then
	1688	write(,) NF_STRERROR(ierr)
	1689	stop 'v'
	1690	endif
	1691
	1692	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
	1693	if(ierr/=NF_NOERR) then
	1694	write(,) NF_STRERROR(ierr)
	1695	stop 'omega'
	1696	endif
	1697
	1698	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
	1699	if(ierr/=NF_NOERR) then
	1700	write(,) NF_STRERROR(ierr)
	1701	stop 'div'
	1702	endif
	1703
	1704	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
	1705	if(ierr/=NF_NOERR) then
	1706	write(,) NF_STRERROR(ierr)
	1707	stop 'T_adv_h'
	1708	endif
	1709
	1710	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
	1711	if(ierr/=NF_NOERR) then
	1712	write(,) NF_STRERROR(ierr)
	1713	stop 'T_adv_v'
	1714	endif
	1715
	1716	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
	1717	if(ierr/=NF_NOERR) then
	1718	write(,) NF_STRERROR(ierr)
	1719	stop 'q_adv_h'
	1720	endif
	1721
	1722	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
	1723	if(ierr/=NF_NOERR) then
	1724	write(,) NF_STRERROR(ierr)
	1725	stop 'q_adv_v'
	1726	endif
	1727
	1728	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
	1729	if(ierr/=NF_NOERR) then
	1730	write(,) NF_STRERROR(ierr)
	1731	stop 's'
	1732	endif
	1733
	1734	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
	1735	if(ierr/=NF_NOERR) then
	1736	write(,) NF_STRERROR(ierr)
	1737	stop 's_adv_h'
	1738	endif
	1739
	1740	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
	1741	if(ierr/=NF_NOERR) then
	1742	write(,) NF_STRERROR(ierr)
	1743	stop 's_adv_v'
	1744	endif
	1745
	1746	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
	1747	if(ierr/=NF_NOERR) then
	1748	write(,) NF_STRERROR(ierr)
	1749	stop 'p_srf_aver'
	1750	endif
	1751
	1752	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
	1753	if(ierr/=NF_NOERR) then
	1754	write(,) NF_STRERROR(ierr)
	1755	stop 'p_srf_center'
	1756	endif
	1757
	1758	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
	1759	if(ierr/=NF_NOERR) then
	1760	write(,) NF_STRERROR(ierr)
	1761	stop 'T_srf'
	1762	endif
	1763
	1764	!dimensions lecture
	1765	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
	1766
	1767	!pressure
	1768	do l=1,ntime
	1769	do k=1,nlevel
	1770	plev(k,l)=lev(k)
	1771	enddo
	1772	enddo
	1773
	1774	#ifdef NC_DOUBLE
	1775	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
	1776	#else
	1777	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
	1778	#endif
	1779	if(ierr/=NF_NOERR) then
	1780	write(,) NF_STRERROR(ierr)
	1781	stop "getvarup"
	1782	endif
	1783	! write(,)'lecture lat ok',lat
	1784
	1785	#ifdef NC_DOUBLE
	1786	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
	1787	#else
	1788	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
	1789	#endif
	1790	if(ierr/=NF_NOERR) then
	1791	write(,) NF_STRERROR(ierr)
	1792	stop "getvarup"
	1793	endif
	1794	! write(,)'lecture lon ok',lon
	1795
	1796	#ifdef NC_DOUBLE
	1797	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
	1798	#else
	1799	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
	1800	#endif
	1801	if(ierr/=NF_NOERR) then
	1802	write(,) NF_STRERROR(ierr)
	1803	stop "getvarup"
	1804	endif
	1805	! write(,)'lecture alt ok',alt
	1806
	1807	#ifdef NC_DOUBLE
	1808	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
	1809	#else
	1810	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
	1811	#endif
	1812	if(ierr/=NF_NOERR) then
	1813	write(,) NF_STRERROR(ierr)
	1814	stop "getvarup"
	1815	endif
	1816	! write(,)'lecture phis ok',phis
	1817
	1818	#ifdef NC_DOUBLE
	1819	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
	1820	#else
	1821	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
	1822	#endif
	1823	if(ierr/=NF_NOERR) then
	1824	write(,) NF_STRERROR(ierr)
	1825	stop "getvarup"
	1826	endif
	1827	! write(,)'lecture T ok'
	1828
	1829	#ifdef NC_DOUBLE
	1830	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
	1831	#else
	1832	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
	1833	#endif
	1834	if(ierr/=NF_NOERR) then
	1835	write(,) NF_STRERROR(ierr)
	1836	stop "getvarup"
	1837	endif
	1838	! write(,)'lecture q ok'
	1839	!q in kg/kg
	1840	do l=1,ntime
	1841	do k=1,nlevel
	1842	q(k,l)=q(k,l)/1000.
	1843	enddo
	1844	enddo
	1845	#ifdef NC_DOUBLE
	1846	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
	1847	#else
	1848	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
	1849	#endif
	1850	if(ierr/=NF_NOERR) then
	1851	write(,) NF_STRERROR(ierr)
	1852	stop "getvarup"
	1853	endif
	1854	! write(,)'lecture u ok'
	1855
	1856	#ifdef NC_DOUBLE
	1857	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
	1858	#else
	1859	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
	1860	#endif
	1861	if(ierr/=NF_NOERR) then
	1862	write(,) NF_STRERROR(ierr)
	1863	stop "getvarup"
	1864	endif
	1865	! write(,)'lecture v ok'
	1866
	1867	#ifdef NC_DOUBLE
	1868	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
	1869	#else
	1870	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
	1871	#endif
	1872	if(ierr/=NF_NOERR) then
	1873	write(,) NF_STRERROR(ierr)
	1874	stop "getvarup"
	1875	endif
	1876	! write(,)'lecture omega ok'
	1877	!omega in mb/hour
	1878	do l=1,ntime
	1879	do k=1,nlevel
	1880	omega(k,l)=omega(k,l)*100./3600.
	1881	enddo
	1882	enddo
	1883
	1884	#ifdef NC_DOUBLE
	1885	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
	1886	#else
	1887	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
	1888	#endif
	1889	if(ierr/=NF_NOERR) then
	1890	write(,) NF_STRERROR(ierr)
	1891	stop "getvarup"
	1892	endif
	1893	! write(,)'lecture div ok'
	1894
	1895	#ifdef NC_DOUBLE
	1896	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
	1897	#else
	1898	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
	1899	#endif
	1900	if(ierr/=NF_NOERR) then
	1901	write(,) NF_STRERROR(ierr)
	1902	stop "getvarup"
	1903	endif
	1904	! write(,)'lecture T_adv_h ok'
	1905	!T adv in K/s
	1906	do l=1,ntime
	1907	do k=1,nlevel
	1908	T_adv_h(k,l)=T_adv_h(k,l)/3600.
	1909	enddo
	1910	enddo
	1911
	1912
	1913	#ifdef NC_DOUBLE
	1914	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
	1915	#else
	1916	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
	1917	#endif
	1918	if(ierr/=NF_NOERR) then
	1919	write(,) NF_STRERROR(ierr)
	1920	stop "getvarup"
	1921	endif
	1922	! write(,)'lecture T_adv_v ok'
	1923	!T adv in K/s
	1924	do l=1,ntime
	1925	do k=1,nlevel
	1926	T_adv_v(k,l)=T_adv_v(k,l)/3600.
	1927	enddo
	1928	enddo
	1929
	1930	#ifdef NC_DOUBLE
	1931	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
	1932	#else
	1933	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
	1934	#endif
	1935	if(ierr/=NF_NOERR) then
	1936	write(,) NF_STRERROR(ierr)
	1937	stop "getvarup"
	1938	endif
	1939	! write(,)'lecture q_adv_h ok'
	1940	!q adv in kg/kg/s
	1941	do l=1,ntime
	1942	do k=1,nlevel
	1943	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
	1944	enddo
	1945	enddo
	1946
	1947
	1948	#ifdef NC_DOUBLE
	1949	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
	1950	#else
	1951	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
	1952	#endif
	1953	if(ierr/=NF_NOERR) then
	1954	write(,) NF_STRERROR(ierr)
	1955	stop "getvarup"
	1956	endif
	1957	! write(,)'lecture q_adv_v ok'
	1958	!q adv in kg/kg/s
	1959	do l=1,ntime
	1960	do k=1,nlevel
	1961	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
	1962	enddo
	1963	enddo
	1964
	1965
	1966	#ifdef NC_DOUBLE
	1967	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
	1968	#else
	1969	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
	1970	#endif
	1971	if(ierr/=NF_NOERR) then
	1972	write(,) NF_STRERROR(ierr)
	1973	stop "getvarup"
	1974	endif
	1975
	1976	#ifdef NC_DOUBLE
	1977	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
	1978	#else
	1979	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
	1980	#endif
	1981	if(ierr/=NF_NOERR) then
	1982	write(,) NF_STRERROR(ierr)
	1983	stop "getvarup"
	1984	endif
	1985
	1986	#ifdef NC_DOUBLE
	1987	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
	1988	#else
	1989	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
	1990	#endif
	1991	if(ierr/=NF_NOERR) then
	1992	write(,) NF_STRERROR(ierr)
	1993	stop "getvarup"
	1994	endif
	1995
	1996	#ifdef NC_DOUBLE
	1997	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
	1998	#else
	1999	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
	2000	#endif
	2001	if(ierr/=NF_NOERR) then
	2002	write(,) NF_STRERROR(ierr)
	2003	stop "getvarup"
	2004	endif
	2005
	2006	#ifdef NC_DOUBLE
	2007	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
	2008	#else
	2009	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
	2010	#endif
	2011	if(ierr/=NF_NOERR) then
	2012	write(,) NF_STRERROR(ierr)
	2013	stop "getvarup"
	2014	endif
	2015
	2016	#ifdef NC_DOUBLE
	2017	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
	2018	#else
	2019	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
	2020	#endif
	2021	if(ierr/=NF_NOERR) then
	2022	write(,) NF_STRERROR(ierr)
	2023	stop "getvarup"
	2024	endif
	2025	! write(,)'lecture T_srf ok', T_srf
	2026
	2027	return
	2028	end subroutine read_twpice
	2029	!=====================================================================
	2030	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
	2031
	2032	! use netcdf
	2033
	2034	implicit none
	2035	#include "netcdf.inc"
	2036	integer nid,ttm,llm
	2037	real*8 time(ttm)
	2038	real*8 lev(llm)
	2039	integer ierr
	2040
	2041	integer timevar,levvar
	2042	integer timelen,levlen
	2043	integer timedimin,levdimin
	2044
	2045	! Control & lecture on dimensions
	2046	! ===============================
	2047	ierr=NF_INQ_DIMID(nid,"time",timedimin)
	2048	ierr=NF_INQ_VARID(nid,"time",timevar)
	2049	if (ierr.NE.NF_NOERR) then
	2050	write(,) 'ERROR: Field <time> is missing'
	2051	stop ""
	2052	endif
	2053	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
	2054
	2055	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
	2056	ierr=NF_INQ_VARID(nid,"lev",levvar)
	2057	if (ierr.NE.NF_NOERR) then
	2058	write(,) 'ERROR: Field <lev> is lacking'
	2059	stop ""
	2060	endif
	2061	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
	2062
	2063	if((timelen/=ttm).or.(levlen/=llm)) then
	2064	write(,) 'ERROR: Not the good lenght for axis'
	2065	write(,) 'longitude: ',timelen,ttm+1
	2066	write(,) 'latitude: ',levlen,llm
	2067	stop ""
	2068	endif
	2069
	2070	!#ifdef NC_DOUBLE
	2071	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
	2072	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
	2073	!#else
	2074	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
	2075	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
	2076	!#endif
	2077
	2078	return
	2079	end
	2080	!=====================================================================
	2081
[2019]	2082	SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof &
	2083	& ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof &
	2084	& ,omega_prof,o3mmr_prof &
	2085	& ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &
	2086	& ,omega_mod,o3mmr_mod,mxcalc)
[2017]	2087
	2088	implicit none
	2089
	2090	#include "dimensions.h"
	2091
[2019]	2092	!-------------------------------------------------------------------------
	2093	! Vertical interpolation of SANDUREF forcing data onto model levels
	2094	!-------------------------------------------------------------------------
[2017]	2095
	2096	integer nlevmax
	2097	parameter (nlevmax=41)
	2098	integer nlev_sandu,mxcalc
	2099	! real play(llm), plev_prof(nlevmax)
	2100	! real t_prof(nlevmax),q_prof(nlevmax)
	2101	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	2102	! real ht_prof(nlevmax),vt_prof(nlevmax)
	2103	! real hq_prof(nlevmax),vq_prof(nlevmax)
	2104
	2105	real play(llm), plev_prof(nlev_sandu)
	2106	real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
	2107	real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
	2108	real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
	2109
	2110	real t_mod(llm),thl_mod(llm),q_mod(llm)
	2111	real u_mod(llm),v_mod(llm), w_mod(llm)
	2112	real omega_mod(llm),o3mmr_mod(llm)
	2113
	2114	integer l,k,k1,k2
	2115	real frac,frac1,frac2,fact
	2116
	2117	do l = 1, llm
	2118
	2119	if (play(l).ge.plev_prof(nlev_sandu)) then
	2120
	2121	mxcalc=l
	2122	k1=0
	2123	k2=0
	2124
	2125	if (play(l).le.plev_prof(1)) then
	2126
	2127	do k = 1, nlev_sandu-1
[2019]	2128	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
[2017]	2129	k1=k
	2130	k2=k+1
	2131	endif
	2132	enddo
	2133
	2134	if (k1.eq.0 .or. k2.eq.0) then
	2135	write(,) 'PB! k1, k2 = ',k1,k2
	2136	write(,) 'l,play(l) = ',l,play(l)/100
	2137	do k = 1, nlev_sandu-1
	2138	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
	2139	enddo
	2140	endif
	2141
	2142	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
	2143	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
	2144	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
	2145	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
	2146	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
	2147	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
	2148	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
[2019]	2149	omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
	2150	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
[2017]	2151
	2152	else !play>plev_prof(1)
	2153
	2154	k1=1
	2155	k2=2
	2156	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
	2157	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
	2158	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
	2159	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
	2160	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
	2161	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
	2162	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
	2163	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
	2164	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
	2165	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
	2166
	2167	endif ! play.le.plev_prof(1)
	2168
	2169	else ! above max altitude of forcing file
	2170
[2019]	2171	!jyg
[2017]	2172	fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
	2173	fact = max(fact,0.) !jyg
	2174	fact = exp(-fact) !jyg
	2175	t_mod(l)= t_prof(nlev_sandu) !jyg
	2176	thl_mod(l)= thl_prof(nlev_sandu) !jyg
	2177	q_mod(l)= q_prof(nlev_sandu)*fact !jyg
	2178	u_mod(l)= u_prof(nlev_sandu)*fact !jyg
	2179	v_mod(l)= v_prof(nlev_sandu)*fact !jyg
	2180	w_mod(l)= w_prof(nlev_sandu)*fact !jyg
	2181	omega_mod(l)= omega_prof(nlev_sandu)*fact !jyg
	2182	o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact !jyg
	2183
	2184	endif ! play
	2185
	2186	enddo ! l
	2187
	2188	do l = 1,llm
	2189	! print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
	2190	! $ l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
	2191	enddo
	2192
	2193	return
	2194	end
	2195	!=====================================================================
[2019]	2196	SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof &
	2197	& ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof &
	2198	& ,w_prof,tke_prof,o3mmr_prof &
	2199	& ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &
	2200	& ,tke_mod,o3mmr_mod,mxcalc)
[2017]	2201
	2202	implicit none
	2203
	2204	#include "dimensions.h"
	2205
[2019]	2206	!-------------------------------------------------------------------------
	2207	! Vertical interpolation of Astex forcing data onto model levels
	2208	!-------------------------------------------------------------------------
[2017]	2209
	2210	integer nlevmax
	2211	parameter (nlevmax=41)
	2212	integer nlev_astex,mxcalc
	2213	! real play(llm), plev_prof(nlevmax)
	2214	! real t_prof(nlevmax),qv_prof(nlevmax)
	2215	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	2216	! real ht_prof(nlevmax),vt_prof(nlevmax)
	2217	! real hq_prof(nlevmax),vq_prof(nlevmax)
	2218
	2219	real play(llm), plev_prof(nlev_astex)
	2220	real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
	2221	real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
	2222	real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
	2223	real qt_prof(nlev_astex),tke_prof(nlev_astex)
	2224
	2225	real t_mod(llm),thl_mod(llm),qv_mod(llm)
	2226	real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
	2227	real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
	2228
	2229	integer l,k,k1,k2
	2230	real frac,frac1,frac2,fact
	2231
	2232	do l = 1, llm
	2233
	2234	if (play(l).ge.plev_prof(nlev_astex)) then
	2235
	2236	mxcalc=l
	2237	k1=0
	2238	k2=0
	2239
	2240	if (play(l).le.plev_prof(1)) then
	2241
	2242	do k = 1, nlev_astex-1
[2019]	2243	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
[2017]	2244	k1=k
	2245	k2=k+1
	2246	endif
	2247	enddo
	2248
	2249	if (k1.eq.0 .or. k2.eq.0) then
	2250	write(,) 'PB! k1, k2 = ',k1,k2
	2251	write(,) 'l,play(l) = ',l,play(l)/100
	2252	do k = 1, nlev_astex-1
	2253	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
	2254	enddo
	2255	endif
	2256
	2257	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
	2258	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
	2259	thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
	2260	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
	2261	ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
	2262	qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
	2263	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
	2264	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
	2265	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
	2266	tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
[2019]	2267	o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
[2017]	2268
	2269	else !play>plev_prof(1)
	2270
	2271	k1=1
	2272	k2=2
	2273	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
	2274	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
	2275	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
	2276	thl_mod(l)= frac1thl_prof(k1) - frac2thl_prof(k2)
	2277	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
	2278	ql_mod(l)= frac1ql_prof(k1) - frac2ql_prof(k2)
	2279	qt_mod(l)= frac1qt_prof(k1) - frac2qt_prof(k2)
	2280	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
	2281	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
	2282	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
	2283	tke_mod(l)= frac1tke_prof(k1) - frac2tke_prof(k2)
	2284	o3mmr_mod(l)= frac1o3mmr_prof(k1) - frac2o3mmr_prof(k2)
	2285
	2286	endif ! play.le.plev_prof(1)
	2287
	2288	else ! above max altitude of forcing file
	2289
[2019]	2290	!jyg
[2017]	2291	fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
	2292	fact = max(fact,0.) !jyg
	2293	fact = exp(-fact) !jyg
	2294	t_mod(l)= t_prof(nlev_astex) !jyg
	2295	thl_mod(l)= thl_prof(nlev_astex) !jyg
	2296	qv_mod(l)= qv_prof(nlev_astex)*fact !jyg
	2297	ql_mod(l)= ql_prof(nlev_astex)*fact !jyg
	2298	qt_mod(l)= qt_prof(nlev_astex)*fact !jyg
	2299	u_mod(l)= u_prof(nlev_astex)*fact !jyg
	2300	v_mod(l)= v_prof(nlev_astex)*fact !jyg
	2301	w_mod(l)= w_prof(nlev_astex)*fact !jyg
	2302	tke_mod(l)= tke_prof(nlev_astex)*fact !jyg
	2303	o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact !jyg
	2304
	2305	endif ! play
	2306
	2307	enddo ! l
	2308
	2309	do l = 1,llm
	2310	! print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
	2311	! $ l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
	2312	enddo
	2313
	2314	return
	2315	end
	2316
	2317	!======================================================================
[2019]	2318	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play &
	2319	& ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &
	2320	& ,dth_dyn,dqh_dyn)
[2017]	2321	implicit none
	2322
[2019]	2323	!-------------------------------------------------------------------------
	2324	! Read RICO forcing data
	2325	!-------------------------------------------------------------------------
[2017]	2326	#include "dimensions.h"
	2327
	2328
	2329	integer nlev_rico
	2330	real ts_rico,ps_rico
	2331	real t_rico(llm),q_rico(llm)
	2332	real u_rico(llm),v_rico(llm)
	2333	real w_rico(llm)
	2334	real dth_dyn(llm)
	2335	real dqh_dyn(llm)
	2336
	2337
	2338	real play(llm),zlay(llm)
	2339
	2340
	2341	real prico(nlev_rico),zrico(nlev_rico)
	2342
	2343	character*80 fich_rico
	2344
	2345	integer k,l
	2346
	2347
	2348	print*,fich_rico
	2349	open(21,file=trim(fich_rico),form='formatted')
	2350	do k=1,llm
	2351	zlay(k)=0.
	2352	enddo
	2353
	2354	read(21,*) ps_rico,ts_rico
	2355	prico(1)=ps_rico
	2356	zrico(1)=0.0
	2357	do l=2,nlev_rico
	2358	read(21,*) k,prico(l),zrico(l)
	2359	enddo
	2360	close(21)
	2361
	2362	do k=1,llm
	2363	do l=1,80
	2364	if(prico(l)>play(k)) then
	2365	if(play(k)>prico(l+1)) then
[2019]	2366	zlay(k)=zrico(l)+(play(k)-prico(l)) * &
[2017]	2367	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
[2019]	2368	else
	2369	zlay(k)=zrico(l)+(play(k)-prico(80))* &
[2017]	2370	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
	2371	endif
	2372	endif
	2373	enddo
	2374	print*,k,zlay(k)
	2375	! U
	2376	if(0 < zlay(k) .and. zlay(k) < 4000) then
	2377	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
	2378	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
[2019]	2379	u_rico(k)= -1.9 + (30.0 + 1.9) / &
	2380	& (12000 - 4000) * (zlay(k) - 4000)
[2017]	2381	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
	2382	u_rico(k)=30.0
	2383	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
[2019]	2384	u_rico(k)=30.0 - (30.0) / &
	2385	& (20000 - 13000) * (zlay(k) - 13000)
[2017]	2386	else
	2387	u_rico(k)=0.0
	2388	endif
	2389
	2390	!Q_v
	2391	if(0 < zlay(k) .and. zlay(k) < 740) then
	2392	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
	2393	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
[2019]	2394	q_rico(k)=13.8 + (2.4 - 13.8) / &
	2395	& (3260 - 740) * (zlay(k) - 740)
[2017]	2396	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
[2019]	2397	q_rico(k)=2.4 + (1.8 - 2.4) / &
	2398	& (4000 - 3260) * (zlay(k) - 3260)
[2017]	2399	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
[2019]	2400	q_rico(k)=1.8 + (0 - 1.8) / &
[2096]	2401	& (9000 - 4000) * (zlay(k) - 4000)
[2017]	2402	else
	2403	q_rico(k)=0.0
	2404	endif
	2405
	2406	!T
	2407	if(0 < zlay(k) .and. zlay(k) < 740) then
	2408	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
	2409	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
[2019]	2410	t_rico(k)=292.0 + (278.0 - 292.0) / &
	2411	& (4000 - 740) * (zlay(k) - 740)
[2017]	2412	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
[2019]	2413	t_rico(k)=278.0 + (203.0 - 278.0) / &
	2414	& (15000 - 4000) * (zlay(k) - 4000)
[2017]	2415	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
[2019]	2416	t_rico(k)=203.0 + (194.0 - 203.0) / &
	2417	& (17500 - 15000)* (zlay(k) - 15000)
[2017]	2418	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
[2019]	2419	t_rico(k)=194.0 + (206.0 - 194.0) / &
	2420	& (20000 - 17500)* (zlay(k) - 17500)
[2017]	2421	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
[2019]	2422	t_rico(k)=206.0 + (270.0 - 206.0) / &
	2423	& (60000 - 20000)* (zlay(k) - 20000)
[2017]	2424	endif
	2425
	2426	! W
	2427	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
	2428	w_rico(k)=- (0.005/2260) * zlay(k)
	2429	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
	2430	w_rico(k)=- 0.005
	2431	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
	2432	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
	2433	else
	2434	w_rico(k)=0.0
	2435	endif
	2436
	2437	! dThrz+dTsw0+dTlw0
	2438	if(0 < zlay(k) .and. zlay(k) < 4000) then
[2019]	2439	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ &
	2440	& (864004000) zlay(k)
[2017]	2441	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
[2019]	2442	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) / &
	2443	& (86400(5000 - 4000)) (zlay(k) - 4000)
[2017]	2444	else
	2445	dth_dyn(k)=0.0
	2446	endif
	2447	! dQhrz
	2448	if(0 < zlay(k) .and. zlay(k) < 3000) then
[2019]	2449	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/ &
	2450	& (864003000) (zlay(k))
[2017]	2451	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
	2452	dqh_dyn(k)=0.345 / 86400
	2453	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
[2019]	2454	dqh_dyn(k)=0.345 / 86400 + &
	2455	& (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
[2017]	2456	else
	2457	dqh_dyn(k)=0.0
	2458	endif
	2459
	2460	!? if(play(k)>6e4) then
	2461	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
	2462	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
	2463	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
	2464	!? *(6e4-play(k))/(6e4-3e4)
	2465	!? else
	2466	!? ratqs0(1,k)=ratqshaut
	2467	!? endif
	2468
	2469	enddo
	2470
	2471	do k=1,llm
	2472	q_rico(k)=q_rico(k)/1e3
	2473	dqh_dyn(k)=dqh_dyn(k)/1e3
	2474	v_rico(k)=-3.8
	2475	enddo
	2476
	2477	return
	2478	end
	2479
	2480	!======================================================================
[2019]	2481	SUBROUTINE interp_sandu_time(day,day1,annee_ref &
	2482	& ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu &
	2483	& ,nlev_sandu,ts_sandu,ts_prof)
[2017]	2484	implicit none
	2485
	2486	!---------------------------------------------------------------------------------------
	2487	! Time interpolation of a 2D field to the timestep corresponding to day
	2488	!
	2489	! day: current julian day (e.g. 717538.2)
	2490	! day1: first day of the simulation
	2491	! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
	2492	! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
	2493	!---------------------------------------------------------------------------------------
	2494	! inputs:
	2495	integer annee_ref
	2496	integer nt_sandu,nlev_sandu
	2497	integer year_ini_sandu
	2498	real day, day1,day_ini_sandu,dt_sandu
	2499	real ts_sandu(nt_sandu)
	2500	! outputs:
	2501	real ts_prof
	2502	! local:
	2503	integer it_sandu1, it_sandu2
	2504	real timeit,time_sandu1,time_sandu2,frac
	2505	! Check that initial day of the simulation consistent with SANDU period:
	2506	if (annee_ref.ne.2006 ) then
	2507	print*,'Pour SANDUREF, annee_ref doit etre 2006 '
	2508	print*,'Changer annee_ref dans run.def'
	2509	stop
	2510	endif
	2511	! if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
	2512	! print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
	2513	! print*,'Changer dayref dans run.def'
	2514	! stop
	2515	! endif
	2516
	2517	! Determine timestep relative to the 1st day of TOGA-COARE:
	2518	! timeit=(day-day1)*86400.
	2519	! if (annee_ref.eq.1992) then
	2520	! timeit=(day-day_ini_sandu)*86400.
	2521	! else
	2522	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
	2523	! endif
	2524	timeit=(day-day_ini_sandu)*86400
	2525
	2526	! Determine the closest observation times:
	2527	it_sandu1=INT(timeit/dt_sandu)+1
	2528	it_sandu2=it_sandu1 + 1
	2529	time_sandu1=(it_sandu1-1)*dt_sandu
	2530	time_sandu2=(it_sandu2-1)*dt_sandu
	2531	print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
[2019]	2532	print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2', &
	2533	& it_sandu1,it_sandu2,time_sandu1,time_sandu2
[2017]	2534
	2535	if (it_sandu1 .ge. nt_sandu) then
[2019]	2536	write(,) 'PB-stop: day, it_sandu1, it_sandu2, timeit: ' &
	2537	& ,day,it_sandu1,it_sandu2,timeit/86400.
[2017]	2538	stop
	2539	endif
	2540
	2541	! time interpolation:
	2542	frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
	2543	frac=max(frac,0.0)
	2544
[2019]	2545	ts_prof = ts_sandu(it_sandu2) &
	2546	& -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
[2017]	2547
[2019]	2548	print*, &
	2549	&'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:', &
	2550	&day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1, &
	2551	&it_sandu2,ts_prof
[2017]	2552
	2553	return
	2554	END
	2555	!=====================================================================
[2019]	2556	!-------------------------------------------------------------------------
	2557	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu, &
	2558	& sens,flat,adv_theta,rad_theta,adv_qt)
[2017]	2559	implicit none
	2560
[2019]	2561	!-------------------------------------------------------------------------
	2562	! Read ARM_CU case forcing data
	2563	!-------------------------------------------------------------------------
[2017]	2564
	2565	integer nlev_armcu,nt_armcu
	2566	real sens(nt_armcu),flat(nt_armcu)
	2567	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
	2568	character*80 fich_armcu
	2569
	2570	integer ip
	2571
	2572	integer iy,im,id,ih,in
	2573
	2574	print*,'nlev_armcu',nlev_armcu
	2575
	2576	open(21,file=trim(fich_armcu),form='formatted')
	2577	read(21,'(a)')
	2578	do ip = 1, nt_armcu
	2579	read(21,'(a)')
	2580	read(21,'(a)')
[2019]	2581	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip), &
	2582	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
	2583	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip), &
	2584	& adv_theta(ip),rad_theta(ip),adv_qt(ip)
[2017]	2585	enddo
	2586	close(21)
	2587
	2588	223 format(5i3,5f8.3)
	2589
	2590	return
	2591	end
	2592
	2593	!=====================================================================
[2019]	2594	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof &
	2595	& ,t_prof,q_prof,u_prof,v_prof,w_prof &
	2596	& ,ht_prof,vt_prof,hq_prof,vq_prof &
	2597	& ,t_mod,q_mod,u_mod,v_mod,w_mod &
	2598	& ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
[2017]	2599
	2600	implicit none
	2601
	2602	#include "dimensions.h"
	2603
[2019]	2604	!-------------------------------------------------------------------------
	2605	! Vertical interpolation of TOGA-COARE forcing data onto model levels
	2606	!-------------------------------------------------------------------------
[2017]	2607
	2608	integer nlevmax
	2609	parameter (nlevmax=41)
	2610	integer nlev_toga,mxcalc
	2611	! real play(llm), plev_prof(nlevmax)
	2612	! real t_prof(nlevmax),q_prof(nlevmax)
	2613	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	2614	! real ht_prof(nlevmax),vt_prof(nlevmax)
	2615	! real hq_prof(nlevmax),vq_prof(nlevmax)
	2616
	2617	real play(llm), plev_prof(nlev_toga)
	2618	real t_prof(nlev_toga),q_prof(nlev_toga)
	2619	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
	2620	real ht_prof(nlev_toga),vt_prof(nlev_toga)
	2621	real hq_prof(nlev_toga),vq_prof(nlev_toga)
	2622
	2623	real t_mod(llm),q_mod(llm)
	2624	real u_mod(llm),v_mod(llm), w_mod(llm)
	2625	real ht_mod(llm),vt_mod(llm)
	2626	real hq_mod(llm),vq_mod(llm)
	2627
	2628	integer l,k,k1,k2
	2629	real frac,frac1,frac2,fact
	2630
	2631	do l = 1, llm
	2632
	2633	if (play(l).ge.plev_prof(nlev_toga)) then
	2634
	2635	mxcalc=l
	2636	k1=0
	2637	k2=0
	2638
	2639	if (play(l).le.plev_prof(1)) then
	2640
	2641	do k = 1, nlev_toga-1
[2019]	2642	if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
[2017]	2643	k1=k
	2644	k2=k+1
	2645	endif
	2646	enddo
	2647
	2648	if (k1.eq.0 .or. k2.eq.0) then
	2649	write(,) 'PB! k1, k2 = ',k1,k2
	2650	write(,) 'l,play(l) = ',l,play(l)/100
	2651	do k = 1, nlev_toga-1
	2652	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
	2653	enddo
	2654	endif
	2655
	2656	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
	2657	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
	2658	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
	2659	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
	2660	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
	2661	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
	2662	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
	2663	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
	2664	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
	2665	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
	2666
	2667	else !play>plev_prof(1)
	2668
	2669	k1=1
	2670	k2=2
	2671	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
	2672	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
	2673	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
	2674	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
	2675	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
	2676	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
	2677	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
	2678	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
	2679	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
	2680	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
	2681	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
	2682
	2683	endif ! play.le.plev_prof(1)
	2684
	2685	else ! above max altitude of forcing file
	2686
[2019]	2687	!jyg
[2017]	2688	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
	2689	fact = max(fact,0.) !jyg
	2690	fact = exp(-fact) !jyg
	2691	t_mod(l)= t_prof(nlev_toga) !jyg
	2692	q_mod(l)= q_prof(nlev_toga)*fact !jyg
	2693	u_mod(l)= u_prof(nlev_toga)*fact !jyg
	2694	v_mod(l)= v_prof(nlev_toga)*fact !jyg
	2695	w_mod(l)= 0.0 !jyg
	2696	ht_mod(l)= ht_prof(nlev_toga) !jyg
	2697	vt_mod(l)= vt_prof(nlev_toga) !jyg
	2698	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
	2699	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
	2700
	2701	endif ! play
	2702
	2703	enddo ! l
	2704
	2705	! do l = 1,llm
	2706	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
	2707	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
	2708	! enddo
	2709
	2710	return
	2711	end
	2712
[2126]	2713	!=====================================================================
[2191]	2714	SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas &
	2715	& ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &
	2716	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
	2717	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
	2718	& ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &
	2719	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
	2720	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
	2721
	2722	implicit none
	2723
	2724	#include "dimensions.h"
	2725
	2726	!-------------------------------------------------------------------------
	2727	! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
	2728	!-------------------------------------------------------------------------
	2729
	2730	integer nlevmax
	2731	parameter (nlevmax=41)
	2732	integer nlev_cas,mxcalc
	2733	! real play(llm), plev_prof(nlevmax)
	2734	! real t_prof(nlevmax),q_prof(nlevmax)
	2735	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	2736	! real ht_prof(nlevmax),vt_prof(nlevmax)
	2737	! real hq_prof(nlevmax),vq_prof(nlevmax)
	2738
	2739	real play(llm), plev_prof_cas(nlev_cas)
	2740	real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
	2741	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
	2742	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
	2743	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
	2744	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
	2745	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
	2746	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
	2747
	2748	real t_mod_cas(llm),q_mod_cas(llm)
	2749	real u_mod_cas(llm),v_mod_cas(llm)
	2750	real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
	2751	real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
	2752	real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
	2753	real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
	2754	real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
	2755
	2756	integer l,k,k1,k2
	2757	real frac,frac1,frac2,fact
	2758
	2759	do l = 1, llm
	2760
	2761	if (play(l).ge.plev_prof_cas(nlev_cas)) then
	2762
	2763	mxcalc=l
	2764	k1=0
	2765	k2=0
	2766
	2767	if (play(l).le.plev_prof_cas(1)) then
	2768
	2769	do k = 1, nlev_cas-1
	2770	if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
	2771	k1=k
	2772	k2=k+1
	2773	endif
	2774	enddo
	2775
	2776	if (k1.eq.0 .or. k2.eq.0) then
	2777	write(,) 'PB! k1, k2 = ',k1,k2
	2778	write(,) 'l,play(l) = ',l,play(l)/100
	2779	do k = 1, nlev_cas-1
	2780	write(,) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
	2781	enddo
	2782	endif
	2783
	2784	frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
	2785	t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
	2786	q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
	2787	u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
	2788	v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
	2789	ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
	2790	vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
	2791	w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
	2792	du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
	2793	hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
	2794	vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
	2795	dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
	2796	hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
	2797	vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
	2798	dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
	2799	ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
	2800	vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
	2801	dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
	2802	hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
	2803	vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
[2933]	2804	dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
[2191]	2805
	2806	else !play>plev_prof_cas(1)
	2807
	2808	k1=1
	2809	k2=2
	2810	frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
	2811	frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
	2812	t_mod_cas(l)= frac1t_prof_cas(k1) - frac2t_prof_cas(k2)
	2813	q_mod_cas(l)= frac1q_prof_cas(k1) - frac2q_prof_cas(k2)
	2814	u_mod_cas(l)= frac1u_prof_cas(k1) - frac2u_prof_cas(k2)
	2815	v_mod_cas(l)= frac1v_prof_cas(k1) - frac2v_prof_cas(k2)
	2816	ug_mod_cas(l)= frac1ug_prof_cas(k1) - frac2ug_prof_cas(k2)
	2817	vg_mod_cas(l)= frac1vg_prof_cas(k1) - frac2vg_prof_cas(k2)
	2818	w_mod_cas(l)= frac1vitw_prof_cas(k1) - frac2vitw_prof_cas(k2)
	2819	du_mod_cas(l)= frac1du_prof_cas(k1) - frac2du_prof_cas(k2)
	2820	hu_mod_cas(l)= frac1hu_prof_cas(k1) - frac2hu_prof_cas(k2)
	2821	vu_mod_cas(l)= frac1vu_prof_cas(k1) - frac2vu_prof_cas(k2)
	2822	dv_mod_cas(l)= frac1dv_prof_cas(k1) - frac2dv_prof_cas(k2)
	2823	hv_mod_cas(l)= frac1hv_prof_cas(k1) - frac2hv_prof_cas(k2)
	2824	vv_mod_cas(l)= frac1vv_prof_cas(k1) - frac2vv_prof_cas(k2)
	2825	dt_mod_cas(l)= frac1dt_prof_cas(k1) - frac2dt_prof_cas(k2)
	2826	ht_mod_cas(l)= frac1ht_prof_cas(k1) - frac2ht_prof_cas(k2)
	2827	vt_mod_cas(l)= frac1vt_prof_cas(k1) - frac2vt_prof_cas(k2)
	2828	dq_mod_cas(l)= frac1dq_prof_cas(k1) - frac2dq_prof_cas(k2)
	2829	hq_mod_cas(l)= frac1hq_prof_cas(k1) - frac2hq_prof_cas(k2)
	2830	vq_mod_cas(l)= frac1vq_prof_cas(k1) - frac2vq_prof_cas(k2)
[2933]	2831	dtrad_mod_cas(l)= frac1dtrad_prof_cas(k1) - frac2dtrad_prof_cas(k2)
[2191]	2832
	2833	endif ! play.le.plev_prof_cas(1)
	2834
	2835	else ! above max altitude of forcing file
	2836
	2837	!jyg
	2838	fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
	2839	fact = max(fact,0.) !jyg
	2840	fact = exp(-fact) !jyg
	2841	t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg
	2842	q_mod_cas(l)= q_prof_cas(nlev_cas)*fact !jyg
	2843	u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg
	2844	v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg
	2845	ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg
	2846	vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg
	2847	w_mod_cas(l)= 0.0 !jyg
	2848	du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
	2849	hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg
	2850	vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg
	2851	dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
	2852	hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg
	2853	vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg
	2854	dt_mod_cas(l)= dt_prof_cas(nlev_cas)
	2855	ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg
	2856	vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg
	2857	dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
	2858	hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg
	2859	vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg
[2933]	2860	dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg
[2191]	2861
	2862	endif ! play
	2863
	2864	enddo ! l
	2865
	2866	! do l = 1,llm
	2867	! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
	2868	! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
	2869	! enddo
	2870
	2871	return
	2872	end
	2873	!*****************************************************************************
	2874	!=====================================================================
[2126]	2875	SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof &
	2876	& ,th_prof,qv_prof,u_prof,v_prof,o3_prof &
	2877	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof &
	2878	& ,th_mod,qv_mod,u_mod,v_mod,o3_mod &
	2879	& ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
	2880
	2881	implicit none
	2882
	2883	#include "dimensions.h"
	2884
	2885	!-------------------------------------------------------------------------
	2886	! Vertical interpolation of Dice forcing data onto model levels
	2887	!-------------------------------------------------------------------------
	2888
	2889	integer nlevmax
	2890	parameter (nlevmax=41)
	2891	integer nlev_dice,mxcalc,nt_dice
	2892
	2893	real play(llm), plev_prof(nlev_dice)
	2894	real th_prof(nlev_dice),qv_prof(nlev_dice)
	2895	real u_prof(nlev_dice),v_prof(nlev_dice)
	2896	real o3_prof(nlev_dice)
	2897	real ht_prof(nlev_dice),hq_prof(nlev_dice)
	2898	real hu_prof(nlev_dice),hv_prof(nlev_dice)
	2899	real w_prof(nlev_dice),omega_prof(nlev_dice)
	2900
	2901	real th_mod(llm),qv_mod(llm)
	2902	real u_mod(llm),v_mod(llm), o3_mod(llm)
	2903	real ht_mod(llm),hq_mod(llm)
	2904	real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
	2905
	2906	integer l,k,k1,k2,kp
	2907	real aa,frac,frac1,frac2,fact
	2908
	2909	do l = 1, llm
	2910
	2911	if (play(l).ge.plev_prof(nlev_dice)) then
	2912
	2913	mxcalc=l
	2914	k1=0
	2915	k2=0
	2916
	2917	if (play(l).le.plev_prof(1)) then
	2918
	2919	do k = 1, nlev_dice-1
	2920	if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
	2921	k1=k
	2922	k2=k+1
	2923	endif
	2924	enddo
	2925
	2926	if (k1.eq.0 .or. k2.eq.0) then
	2927	write(,) 'PB! k1, k2 = ',k1,k2
	2928	write(,) 'l,play(l) = ',l,play(l)/100
	2929	do k = 1, nlev_dice-1
	2930	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
	2931	enddo
	2932	endif
	2933
	2934	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
	2935	th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
	2936	qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
	2937	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
	2938	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
	2939	o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
	2940	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
	2941	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
	2942	hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
	2943	hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
	2944	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
	2945	omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
	2946
	2947	else !play>plev_prof(1)
	2948
	2949	k1=1
	2950	k2=2
	2951	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
	2952	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
	2953	th_mod(l)= frac1th_prof(k1) - frac2th_prof(k2)
	2954	qv_mod(l)= frac1qv_prof(k1) - frac2qv_prof(k2)
	2955	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
	2956	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
	2957	o3_mod(l)= frac1o3_prof(k1) - frac2o3_prof(k2)
	2958	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
	2959	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
	2960	hu_mod(l)= frac1hu_prof(k1) - frac2hu_prof(k2)
	2961	hv_mod(l)= frac1hv_prof(k1) - frac2hv_prof(k2)
	2962	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
	2963	omega_mod(l)= frac1omega_prof(k1) - frac2omega_prof(k2)
	2964
	2965	endif ! play.le.plev_prof(1)
	2966
	2967	else ! above max altitude of forcing file
	2968
	2969	!jyg
	2970	fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
	2971	fact = max(fact,0.) !jyg
	2972	fact = exp(-fact) !jyg
	2973	th_mod(l)= th_prof(nlev_dice) !jyg
	2974	qv_mod(l)= qv_prof(nlev_dice)*fact !jyg
	2975	u_mod(l)= u_prof(nlev_dice)*fact !jyg
	2976	v_mod(l)= v_prof(nlev_dice)*fact !jyg
	2977	o3_mod(l)= o3_prof(nlev_dice)*fact !jyg
	2978	ht_mod(l)= ht_prof(nlev_dice) !jyg
	2979	hq_mod(l)= hq_prof(nlev_dice)*fact !jyg
	2980	hu_mod(l)= hu_prof(nlev_dice) !jyg
	2981	hv_mod(l)= hv_prof(nlev_dice) !jyg
	2982	w_mod(l)= 0. !jyg
	2983	omega_mod(l)= 0. !jyg
	2984
	2985	endif ! play
	2986
	2987	enddo ! l
	2988
	2989	! do l = 1,llm
	2990	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
	2991	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
	2992	! enddo
	2993
	2994	return
	2995	end
	2996
[2017]	2997	!======================================================================
[2019]	2998	SUBROUTINE interp_astex_time(day,day1,annee_ref &
	2999	& ,year_ini_astex,day_ini_astex,nt_astex,dt_astex &
	3000	& ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &
	3001	& ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &
	3002	& ,ufa_prof,vfa_prof)
[2017]	3003	implicit none
	3004
	3005	!---------------------------------------------------------------------------------------
	3006	! Time interpolation of a 2D field to the timestep corresponding to day
	3007	!
	3008	! day: current julian day (e.g. 717538.2)
	3009	! day1: first day of the simulation
	3010	! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
	3011	! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
	3012	!---------------------------------------------------------------------------------------
	3013
	3014	! inputs:
	3015	integer annee_ref
	3016	integer nt_astex,nlev_astex
	3017	integer year_ini_astex
	3018	real day, day1,day_ini_astex,dt_astex
	3019	real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
	3020	real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
	3021	! outputs:
	3022	real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
	3023	! local:
	3024	integer it_astex1, it_astex2
	3025	real timeit,time_astex1,time_astex2,frac
	3026
	3027	! Check that initial day of the simulation consistent with ASTEX period:
	3028	if (annee_ref.ne.1992 ) then
	3029	print*,'Pour Astex, annee_ref doit etre 1992 '
	3030	print*,'Changer annee_ref dans run.def'
	3031	stop
	3032	endif
	3033	if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
	3034	print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
	3035	print*,'Changer dayref dans run.def'
	3036	stop
	3037	endif
	3038
	3039	! Determine timestep relative to the 1st day of TOGA-COARE:
	3040	! timeit=(day-day1)*86400.
	3041	! if (annee_ref.eq.1992) then
	3042	! timeit=(day-day_ini_astex)*86400.
	3043	! else
	3044	! timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
	3045	! endif
	3046	timeit=(day-day_ini_astex)*86400
	3047
	3048	! Determine the closest observation times:
	3049	it_astex1=INT(timeit/dt_astex)+1
	3050	it_astex2=it_astex1 + 1
	3051	time_astex1=(it_astex1-1)*dt_astex
	3052	time_astex2=(it_astex2-1)*dt_astex
	3053	print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
[2019]	3054	print *,'it_astex1,it_astex2,time_astex1,time_astex2', &
	3055	& it_astex1,it_astex2,time_astex1,time_astex2
[2017]	3056
	3057	if (it_astex1 .ge. nt_astex) then
[2019]	3058	write(,) 'PB-stop: day, it_astex1, it_astex2, timeit: ' &
	3059	& ,day,it_astex1,it_astex2,timeit/86400.
[2017]	3060	stop
	3061	endif
	3062
	3063	! time interpolation:
	3064	frac=(time_astex2-timeit)/(time_astex2-time_astex1)
	3065	frac=max(frac,0.0)
	3066
[2019]	3067	div_prof = div_astex(it_astex2) &
	3068	& -frac*(div_astex(it_astex2)-div_astex(it_astex1))
	3069	ts_prof = ts_astex(it_astex2) &
	3070	& -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
	3071	ug_prof = ug_astex(it_astex2) &
	3072	& -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
	3073	vg_prof = vg_astex(it_astex2) &
	3074	& -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
	3075	ufa_prof = ufa_astex(it_astex2) &
	3076	& -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
	3077	vfa_prof = vfa_astex(it_astex2) &
	3078	& -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
[2017]	3079
[2019]	3080	print*, &
	3081	&'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:', &
	3082	&day,annee_ref,day_ini_astex,timeit/86400.,it_astex1, &
	3083	&it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
[2017]	3084
	3085	return
	3086	END
	3087
	3088	!======================================================================
[2019]	3089	SUBROUTINE interp_toga_time(day,day1,annee_ref &
	3090	& ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga &
	3091	& ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &
	3092	& ,ht_toga,vt_toga,hq_toga,vq_toga &
	3093	& ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &
	3094	& ,ht_prof,vt_prof,hq_prof,vq_prof)
[2017]	3095	implicit none
	3096
	3097	!---------------------------------------------------------------------------------------
	3098	! Time interpolation of a 2D field to the timestep corresponding to day
	3099	!
	3100	! day: current julian day (e.g. 717538.2)
	3101	! day1: first day of the simulation
	3102	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
	3103	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
	3104	!---------------------------------------------------------------------------------------
	3105
	3106	#include "compar1d.h"
	3107
	3108	! inputs:
	3109	integer annee_ref
	3110	integer nt_toga,nlev_toga
	3111	integer year_ini_toga
	3112	real day, day1,day_ini_toga,dt_toga
	3113	real ts_toga(nt_toga)
	3114	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
	3115	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
	3116	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
	3117	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
	3118	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
	3119	! outputs:
	3120	real ts_prof
	3121	real plev_prof(nlev_toga),t_prof(nlev_toga)
	3122	real q_prof(nlev_toga),u_prof(nlev_toga)
	3123	real v_prof(nlev_toga),w_prof(nlev_toga)
	3124	real ht_prof(nlev_toga),vt_prof(nlev_toga)
	3125	real hq_prof(nlev_toga),vq_prof(nlev_toga)
	3126	! local:
	3127	integer it_toga1, it_toga2,k
	3128	real timeit,time_toga1,time_toga2,frac
	3129
	3130
	3131	if (forcing_type.eq.2) then
	3132	! Check that initial day of the simulation consistent with TOGA-COARE period:
	3133	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
	3134	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
	3135	print*,'Changer annee_ref dans run.def'
	3136	stop
	3137	endif
	3138	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
[2716]	3139	print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'
[2017]	3140	print*,'Changer dayref dans run.def'
	3141	stop
	3142	endif
	3143	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
	3144	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
	3145	print*,'Changer dayref ou nday dans run.def'
	3146	stop
	3147	endif
	3148
	3149	else if (forcing_type.eq.4) then
	3150
	3151	! Check that initial day of the simulation consistent with TWP-ICE period:
	3152	if (annee_ref.ne.2006) then
	3153	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
	3154	print*,'Changer annee_ref dans run.def'
	3155	stop
	3156	endif
	3157	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
	3158	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
	3159	print*,'Changer dayref dans run.def'
	3160	stop
	3161	endif
	3162	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
	3163	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
	3164	print*,'Changer dayref ou nday dans run.def'
	3165	stop
	3166	endif
	3167
	3168	endif
	3169
	3170	! Determine timestep relative to the 1st day of TOGA-COARE:
	3171	! timeit=(day-day1)*86400.
	3172	! if (annee_ref.eq.1992) then
	3173	! timeit=(day-day_ini_toga)*86400.
	3174	! else
	3175	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
	3176	! endif
	3177	timeit=(day-day_ini_toga)*86400
	3178
	3179	! Determine the closest observation times:
	3180	it_toga1=INT(timeit/dt_toga)+1
	3181	it_toga2=it_toga1 + 1
	3182	time_toga1=(it_toga1-1)*dt_toga
	3183	time_toga2=(it_toga2-1)*dt_toga
	3184
	3185	if (it_toga1 .ge. nt_toga) then
[2019]	3186	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: ' &
	3187	& ,day,it_toga1,it_toga2,timeit/86400.
[2017]	3188	stop
	3189	endif
	3190
	3191	! time interpolation:
	3192	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
	3193	frac=max(frac,0.0)
	3194
[2019]	3195	ts_prof = ts_toga(it_toga2) &
	3196	& -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
[2017]	3197
	3198	! print*,
	3199	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
	3200	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
	3201
	3202	do k=1,nlev_toga
[2019]	3203	plev_prof(k) = 100.*(plev_toga(k,it_toga2) &
	3204	& -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
	3205	t_prof(k) = t_toga(k,it_toga2) &
	3206	& -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
	3207	q_prof(k) = q_toga(k,it_toga2) &
	3208	& -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
	3209	u_prof(k) = u_toga(k,it_toga2) &
	3210	& -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
	3211	v_prof(k) = v_toga(k,it_toga2) &
	3212	& -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
	3213	w_prof(k) = w_toga(k,it_toga2) &
	3214	& -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
	3215	ht_prof(k) = ht_toga(k,it_toga2) &
	3216	& -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
	3217	vt_prof(k) = vt_toga(k,it_toga2) &
	3218	& -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
	3219	hq_prof(k) = hq_toga(k,it_toga2) &
	3220	& -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
	3221	vq_prof(k) = vq_toga(k,it_toga2) &
	3222	& -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
[2017]	3223	enddo
	3224
	3225	return
	3226	END
	3227
	3228	!======================================================================
[2126]	3229	SUBROUTINE interp_dice_time(day,day1,annee_ref &
	3230	& ,year_ini_dice,day_ini_dice,nt_dice,dt_dice &
	3231	& ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice &
	3232	& ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice &
	3233	& ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
	3234	& ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof &
	3235	& ,ustar_prof,psurf_prof,ug_prof,vg_prof &
	3236	& ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
	3237	implicit none
	3238
	3239	!---------------------------------------------------------------------------------------
	3240	! Time interpolation of a 2D field to the timestep corresponding to day
	3241	!
	3242	! day: current julian day (e.g. 717538.2)
	3243	! day1: first day of the simulation
	3244	! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
	3245	! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
	3246	!---------------------------------------------------------------------------------------
	3247
	3248	#include "compar1d.h"
	3249
	3250	! inputs:
	3251	integer annee_ref
	3252	integer nt_dice,nlev_dice
	3253	integer year_ini_dice
	3254	real day, day1,day_ini_dice,dt_dice
	3255	real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
	3256	real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
	3257	real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
	3258	real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
	3259	real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
	3260	real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
	3261	! outputs:
	3262	real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
	3263	real ustar_prof,psurf_prof,ug_prof,vg_prof
	3264	real ht_prof(nlev_dice),hq_prof(nlev_dice)
	3265	real hu_prof(nlev_dice),hv_prof(nlev_dice)
	3266	real w_prof(nlev_dice),omega_prof(nlev_dice)
	3267	! local:
	3268	integer it_dice1, it_dice2,k
	3269	real timeit,time_dice1,time_dice2,frac
	3270
	3271
	3272	if (forcing_type.eq.7) then
	3273	! Check that initial day of the simulation consistent with Dice period:
	3274	print *,'annee_ref=',annee_ref
	3275	print *,'day1=',day1
	3276	print *,'day_ini_dice=',day_ini_dice
	3277	if (annee_ref.ne.1999) then
	3278	print*,'Pour Dice, annee_ref doit etre 1999'
	3279	print*,'Changer annee_ref dans run.def'
	3280	stop
	3281	endif
	3282	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
	3283	print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
	3284	print*,'Changer dayref dans run.def',day1,day_ini_dice
	3285	stop
	3286	endif
	3287	if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
	3288	print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
	3289	print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
	3290	stop
	3291	endif
	3292
	3293	endif
	3294
	3295	! Determine timestep relative to the 1st day of TOGA-COARE:
	3296	! timeit=(day-day1)*86400.
	3297	! if (annee_ref.eq.1992) then
	3298	! timeit=(day-day_ini_dice)*86400.
	3299	! else
	3300	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
	3301	! endif
	3302	timeit=(day-day_ini_dice)*86400
	3303
	3304	! Determine the closest observation times:
	3305	it_dice1=INT(timeit/dt_dice)+1
	3306	it_dice2=it_dice1 + 1
	3307	time_dice1=(it_dice1-1)*dt_dice
	3308	time_dice2=(it_dice2-1)*dt_dice
	3309
	3310	if (it_dice1 .ge. nt_dice) then
	3311	write(,) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
	3312	stop
	3313	endif
	3314
	3315	! time interpolation:
	3316	frac=(time_dice2-timeit)/(time_dice2-time_dice1)
	3317	frac=max(frac,0.0)
	3318
	3319	shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
	3320	lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
	3321	lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
	3322	swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
	3323	tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
	3324	ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
	3325	psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
	3326	ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
	3327	vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
	3328
	3329	! print*,
	3330	! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
	3331	! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
	3332
	3333	do k=1,nlev_dice
	3334	ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
	3335	hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
	3336	hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
	3337	hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
	3338	w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
	3339	omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
	3340	enddo
	3341
	3342	return
	3343	END
	3344
	3345	!======================================================================
[2672]	3346	SUBROUTINE interp_gabls4_time(day,day1,annee_ref &
	3347	& ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 &
	3348	& ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 &
	3349	& ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
	3350	implicit none
	3351
	3352	!---------------------------------------------------------------------------------------
	3353	! Time interpolation of a 2D field to the timestep corresponding to day
	3354	!
	3355	! day: current julian day
	3356	! day1: first day of the simulation
	3357	! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
	3358	! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
	3359	!---------------------------------------------------------------------------------------
	3360
	3361	#include "compar1d.h"
	3362
	3363	! inputs:
	3364	integer annee_ref
	3365	integer nt_gabls4,nlev_gabls4
	3366	integer year_ini_gabls4
	3367	real day, day1,day_ini_gabls4,dt_gabls4
	3368	real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
	3369	real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
	3370	real tg_gabls4(nt_gabls4), tg_prof
	3371	! outputs:
	3372	real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
	3373	real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
	3374	! local:
	3375	integer it_gabls41, it_gabls42,k
	3376	real timeit,time_gabls41,time_gabls42,frac
	3377
	3378
	3379
	3380	! Check that initial day of the simulation consistent with gabls4 period:
	3381	if (forcing_type.eq.8 ) then
	3382	print *,'annee_ref=',annee_ref
	3383	print *,'day1=',day1
	3384	print *,'day_ini_gabls4=',day_ini_gabls4
	3385	if (annee_ref.ne.2009) then
	3386	print*,'Pour gabls4, annee_ref doit etre 2009'
	3387	print*,'Changer annee_ref dans run.def'
	3388	stop
	3389	endif
	3390	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
	3391	print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
	3392	print*,'Changer dayref dans run.def',day1,day_ini_gabls4
	3393	stop
	3394	endif
	3395	if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
	3396	print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
	3397	print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
	3398	stop
	3399	endif
	3400	endif
	3401
	3402	timeit=(day-day_ini_gabls4)*86400
	3403	print *,'day,day_ini_gabls4=',day,day_ini_gabls4
	3404	print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
	3405
	3406	! Determine the closest observation times:
	3407	it_gabls41=INT(timeit/dt_gabls4)+1
	3408	it_gabls42=it_gabls41 + 1
	3409	time_gabls41=(it_gabls41-1)*dt_gabls4
	3410	time_gabls42=(it_gabls42-1)*dt_gabls4
	3411
	3412	if (it_gabls41 .ge. nt_gabls4) then
	3413	write(,) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
	3414	stop
	3415	endif
	3416
	3417	! time interpolation:
	3418	frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
	3419	frac=max(frac,0.0)
	3420
	3421
	3422	do k=1,nlev_gabls4
	3423	ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
	3424	vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
	3425	ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
	3426	hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
	3427	enddo
	3428	tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
	3429	return
	3430	END
	3431
	3432	!======================================================================
[2019]	3433	SUBROUTINE interp_armcu_time(day,day1,annee_ref &
	3434	& ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu &
	3435	& ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu &
	3436	& ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
[2017]	3437	implicit none
	3438
	3439	!---------------------------------------------------------------------------------------
	3440	! Time interpolation of a 2D field to the timestep corresponding to day
	3441	!
	3442	! day: current julian day (e.g. 717538.2)
	3443	! day1: first day of the simulation
	3444	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
	3445	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
	3446	! fs= sensible flux
	3447	! fl= latent flux
	3448	! at,rt,aqt= advective and radiative tendencies
	3449	!---------------------------------------------------------------------------------------
	3450
	3451	! inputs:
	3452	integer annee_ref
	3453	integer nt_armcu,nlev_armcu
	3454	integer year_ini_armcu
	3455	real day, day1,day_ini_armcu,dt_armcu
	3456	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
	3457	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
	3458	! outputs:
	3459	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
	3460	! local:
	3461	integer it_armcu1, it_armcu2,k
	3462	real timeit,time_armcu1,time_armcu2,frac
	3463
	3464	! Check that initial day of the simulation consistent with ARMCU period:
	3465	if (annee_ref.ne.1997 ) then
	3466	print*,'Pour ARMCU, annee_ref doit etre 1997 '
	3467	print*,'Changer annee_ref dans run.def'
	3468	stop
	3469	endif
	3470
	3471	timeit=(day-day_ini_armcu)*86400
	3472
	3473	! Determine the closest observation times:
	3474	it_armcu1=INT(timeit/dt_armcu)+1
	3475	it_armcu2=it_armcu1 + 1
	3476	time_armcu1=(it_armcu1-1)*dt_armcu
	3477	time_armcu2=(it_armcu2-1)*dt_armcu
	3478	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
[2019]	3479	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', &
	3480	& it_armcu1,it_armcu2,time_armcu1,time_armcu2
[2017]	3481
	3482	if (it_armcu1 .ge. nt_armcu) then
[2019]	3483	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' &
	3484	& ,day,it_armcu1,it_armcu2,timeit/86400.
[2017]	3485	stop
	3486	endif
	3487
	3488	! time interpolation:
	3489	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
	3490	frac=max(frac,0.0)
	3491
[2019]	3492	fs_prof = fs_armcu(it_armcu2) &
	3493	& -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
	3494	fl_prof = fl_armcu(it_armcu2) &
	3495	& -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
	3496	at_prof = at_armcu(it_armcu2) &
	3497	& -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
	3498	rt_prof = rt_armcu(it_armcu2) &
	3499	& -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
	3500	aqt_prof = aqt_armcu(it_armcu2) &
	3501	& -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
[2017]	3502
[2019]	3503	print*, &
	3504	&'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', &
	3505	&day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, &
	3506	&it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
[2017]	3507
	3508	return
	3509	END
	3510
	3511	!=====================================================================
[2019]	3512	subroutine readprofiles(nlev_max,kmax,ntrac,height, &
	3513	& thlprof,qtprof,uprof, &
	3514	& vprof,e12prof,ugprof,vgprof, &
	3515	& wfls,dqtdxls,dqtdyls,dqtdtls, &
	3516	& thlpcar,tracer,nt1,nt2)
[2017]	3517	implicit none
	3518
	3519	integer nlev_max,kmax,kmax2,ntrac
	3520	logical :: llesread = .true.
	3521
[2019]	3522	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), &
	3523	& uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), &
	3524	& ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), &
	3525	& dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), &
	3526	& thlpcar(nlev_max),tracer(nlev_max,ntrac)
[2017]	3527
[2565]	3528	real height1(nlev_max)
	3529
[2017]	3530	integer, parameter :: ilesfile=1
	3531	integer :: ierr,k,itrac,nt1,nt2
	3532
	3533	if(.not.(llesread)) return
	3534
	3535	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
	3536	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	3537	read (ilesfile,*) kmax
	3538	do k=1,kmax
[2565]	3539	read (ilesfile,*) height1(k),thlprof(k),qtprof (k), &
[2019]	3540	& uprof (k),vprof (k),e12prof(k)
[2017]	3541	enddo
	3542	close(ilesfile)
	3543
	3544	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
	3545	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
	3546	read (ilesfile,*) kmax2
	3547	if (kmax .ne. kmax2) then
	3548	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
	3549	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
	3550	stop 'lecture profiles'
	3551	endif
	3552	do k=1,kmax
[2019]	3553	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &
	3554	& dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
[2017]	3555	end do
[2565]	3556	do k=1,kmax
	3557	if (height(k) .ne. height1(k)) then
	3558	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
	3559	print *, 'les niveaux different : ',k,height1(k), height(k)
	3560	stop
	3561	endif
	3562	end do
[2017]	3563	close(ilesfile)
	3564
	3565	open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
	3566	if (ierr /= 0) then
	3567	print*,'WARNING : trac.inp does not exist'
	3568	else
	3569	read (ilesfile,*) kmax2,nt1,nt2
	3570	if (nt2>ntrac) then
	3571	stop'Augmenter le nombre de traceurs dans traceur.def'
	3572	endif
	3573	if (kmax .ne. kmax2) then
	3574	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
	3575	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
	3576	stop 'lecture profiles'
	3577	endif
	3578	do k=1,kmax
	3579	read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
	3580	end do
	3581	close(ilesfile)
	3582	endif
	3583
	3584	return
	3585	end
	3586	!======================================================================
[2019]	3587	subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &
	3588	& thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
[2017]	3589	!======================================================================
	3590	implicit none
	3591
	3592	integer nlev_max,kmax
	3593	logical :: llesread = .true.
	3594
[2019]	3595	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
	3596	real thlprof(nlev_max)
	3597	real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
	3598	real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
[2017]	3599
	3600	integer, parameter :: ilesfile=1
	3601	integer :: k,ierr
	3602
	3603	if(.not.(llesread)) return
	3604
	3605	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
	3606	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	3607	read (ilesfile,*) kmax
	3608	do k=1,kmax
[2019]	3609	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
	3610	& qprof (k),uprof(k), vprof(k), wprof(k), &
	3611	& omega (k),o3mmr(k)
[2017]	3612	enddo
	3613	close(ilesfile)
	3614
	3615	return
	3616	end
	3617
	3618	!======================================================================
[2019]	3619	subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &
	3620	& thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
[2017]	3621	!======================================================================
	3622	implicit none
	3623
	3624	integer nlev_max,kmax
	3625	logical :: llesread = .true.
	3626
[2019]	3627	real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &
	3628	& thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &
	3629	& qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &
	3630	& wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
[2017]	3631
	3632	integer, parameter :: ilesfile=1
	3633	integer :: ierr,k
	3634
	3635	if(.not.(llesread)) return
	3636
	3637	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
	3638	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	3639	read (ilesfile,*) kmax
	3640	do k=1,kmax
[2019]	3641	read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
	3642	& qvprof (k),qlprof (k),qtprof (k), &
	3643	& uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
[2017]	3644	enddo
	3645	close(ilesfile)
	3646
	3647	return
	3648	end
	3649
	3650
	3651
	3652	!======================================================================
[2019]	3653	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &
	3654	& vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
[2017]	3655	!======================================================================
	3656	implicit none
	3657
	3658	integer nlev_max,kmax
	3659	logical :: llesread = .true.
	3660
[2019]	3661	real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
	3662	real thetaprof(nlev_max),rvprof(nlev_max)
	3663	real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
	3664	real aprof(nlev_max+1),bprof(nlev_max+1)
[2017]	3665
	3666	integer, parameter :: ilesfile=1
	3667	integer, parameter :: ifile=2
	3668	integer :: ierr,jtot,k
	3669
	3670	if(.not.(llesread)) return
	3671
	3672	! Read profiles at full levels
	3673	IF(nlev_max.EQ.19) THEN
	3674	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
	3675	print *,'On ouvre prof.inp.19'
	3676	ELSE
	3677	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
	3678	print *,'On ouvre prof.inp.40'
	3679	ENDIF
	3680	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	3681	read (ilesfile,*) kmax
	3682	do k=1,kmax
[2019]	3683	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &
	3684	& thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
[2017]	3685	enddo
	3686	close(ilesfile)
	3687
	3688	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
	3689	IF(nlev_max.EQ.19) THEN
	3690	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
	3691	print *,'On ouvre proh.inp.19'
	3692	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
	3693	ELSE
	3694	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
	3695	print *,'On ouvre proh.inp.40'
	3696	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
	3697	ENDIF
	3698	read (ifile,*) kmax
	3699	do k=1,kmax
	3700	read (ifile,*) jtot,aprof(k),bprof(k)
	3701	enddo
	3702	close(ifile)
	3703
	3704	return
	3705	end
	3706
	3707	!=====================================================================
[2019]	3708	subroutine read_fire(fich_fire,nlevel,ntime &
	3709	& ,zz,thl,qt,u,v,tke &
	3710	& ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
[2017]	3711
	3712	!program reading forcings of the FIRE case study
	3713
	3714
	3715	implicit none
	3716
	3717	#include "netcdf.inc"
	3718
	3719	integer ntime,nlevel
	3720	character*80 :: fich_fire
	3721	real*8 zz(nlevel)
	3722
	3723	real*8 thl(nlevel)
	3724	real*8 qt(nlevel),u(nlevel)
	3725	real*8 v(nlevel),tke(nlevel)
	3726	real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
	3727	real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
	3728	real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
	3729
	3730	integer nid, ierr
	3731	integer nbvar3d
	3732	parameter(nbvar3d=30)
	3733	integer var3didin(nbvar3d)
	3734
	3735	ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
	3736	if (ierr.NE.NF_NOERR) then
	3737	write(,) 'ERROR: Pb opening forcings nc file '
	3738	write(,) NF_STRERROR(ierr)
	3739	stop ""
	3740	endif
	3741
	3742
	3743	ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
	3744	if(ierr/=NF_NOERR) then
	3745	write(,) NF_STRERROR(ierr)
	3746	stop 'lev'
	3747	endif
	3748
	3749
	3750	ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
	3751	if(ierr/=NF_NOERR) then
	3752	write(,) NF_STRERROR(ierr)
	3753	stop 'temp'
	3754	endif
	3755
	3756	ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
	3757	if(ierr/=NF_NOERR) then
	3758	write(,) NF_STRERROR(ierr)
	3759	stop 'qv'
	3760	endif
	3761
	3762	ierr=NF_INQ_VARID(nid,"u",var3didin(4))
	3763	if(ierr/=NF_NOERR) then
	3764	write(,) NF_STRERROR(ierr)
	3765	stop 'u'
	3766	endif
	3767
	3768	ierr=NF_INQ_VARID(nid,"v",var3didin(5))
	3769	if(ierr/=NF_NOERR) then
	3770	write(,) NF_STRERROR(ierr)
	3771	stop 'v'
	3772	endif
	3773
	3774	ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
	3775	if(ierr/=NF_NOERR) then
	3776	write(,) NF_STRERROR(ierr)
	3777	stop 'tke'
	3778	endif
	3779
	3780	ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
	3781	if(ierr/=NF_NOERR) then
	3782	write(,) NF_STRERROR(ierr)
	3783	stop 'ug'
	3784	endif
	3785
	3786	ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
	3787	if(ierr/=NF_NOERR) then
	3788	write(,) NF_STRERROR(ierr)
	3789	stop 'vg'
	3790	endif
	3791
	3792	ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
	3793	if(ierr/=NF_NOERR) then
	3794	write(,) NF_STRERROR(ierr)
	3795	stop 'wls'
	3796	endif
	3797
	3798	ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
	3799	if(ierr/=NF_NOERR) then
	3800	write(,) NF_STRERROR(ierr)
	3801	stop 'dqtdx'
	3802	endif
	3803
	3804	ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
	3805	if(ierr/=NF_NOERR) then
	3806	write(,) NF_STRERROR(ierr)
	3807	stop 'dqtdy'
	3808	endif
	3809
	3810	ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
	3811	if(ierr/=NF_NOERR) then
	3812	write(,) NF_STRERROR(ierr)
	3813	stop 'dqtdt'
	3814	endif
	3815
	3816	ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
	3817	if(ierr/=NF_NOERR) then
	3818	write(,) NF_STRERROR(ierr)
	3819	stop 'thl_rad'
	3820	endif
	3821	!dimensions lecture
	3822	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
	3823
	3824	#ifdef NC_DOUBLE
	3825	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
	3826	#else
	3827	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
	3828	#endif
	3829	if(ierr/=NF_NOERR) then
	3830	write(,) NF_STRERROR(ierr)
	3831	stop "getvarup"
	3832	endif
	3833	! write(,)'lecture z ok',zz
	3834
	3835	#ifdef NC_DOUBLE
	3836	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
	3837	#else
	3838	ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
	3839	#endif
	3840	if(ierr/=NF_NOERR) then
	3841	write(,) NF_STRERROR(ierr)
	3842	stop "getvarup"
	3843	endif
	3844	! write(,)'lecture thl ok',thl
	3845
	3846	#ifdef NC_DOUBLE
	3847	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
	3848	#else
	3849	ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
	3850	#endif
	3851	if(ierr/=NF_NOERR) then
	3852	write(,) NF_STRERROR(ierr)
	3853	stop "getvarup"
	3854	endif
	3855	! write(,)'lecture qt ok',qt
	3856
	3857	#ifdef NC_DOUBLE
	3858	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
	3859	#else
	3860	ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
	3861	#endif
	3862	if(ierr/=NF_NOERR) then
	3863	write(,) NF_STRERROR(ierr)
	3864	stop "getvarup"
	3865	endif
	3866	! write(,)'lecture u ok',u
	3867
	3868	#ifdef NC_DOUBLE
	3869	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
	3870	#else
	3871	ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
	3872	#endif
	3873	if(ierr/=NF_NOERR) then
	3874	write(,) NF_STRERROR(ierr)
	3875	stop "getvarup"
	3876	endif
	3877	! write(,)'lecture v ok',v
	3878
	3879	#ifdef NC_DOUBLE
	3880	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
	3881	#else
	3882	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
	3883	#endif
	3884	if(ierr/=NF_NOERR) then
	3885	write(,) NF_STRERROR(ierr)
	3886	stop "getvarup"
	3887	endif
	3888	! write(,)'lecture tke ok',tke
	3889
	3890	#ifdef NC_DOUBLE
	3891	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
	3892	#else
	3893	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
	3894	#endif
	3895	if(ierr/=NF_NOERR) then
	3896	write(,) NF_STRERROR(ierr)
	3897	stop "getvarup"
	3898	endif
	3899	! write(,)'lecture ug ok',ug
	3900
	3901	#ifdef NC_DOUBLE
	3902	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
	3903	#else
	3904	ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
	3905	#endif
	3906	if(ierr/=NF_NOERR) then
	3907	write(,) NF_STRERROR(ierr)
	3908	stop "getvarup"
	3909	endif
	3910	! write(,)'lecture vg ok',vg
	3911
	3912	#ifdef NC_DOUBLE
	3913	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
	3914	#else
	3915	ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
	3916	#endif
	3917	if(ierr/=NF_NOERR) then
	3918	write(,) NF_STRERROR(ierr)
	3919	stop "getvarup"
	3920	endif
	3921	! write(,)'lecture wls ok',wls
	3922
	3923	#ifdef NC_DOUBLE
	3924	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
	3925	#else
	3926	ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
	3927	#endif
	3928	if(ierr/=NF_NOERR) then
	3929	write(,) NF_STRERROR(ierr)
	3930	stop "getvarup"
	3931	endif
	3932	! write(,)'lecture dqtdx ok',dqtdx
	3933
	3934	#ifdef NC_DOUBLE
	3935	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
	3936	#else
	3937	ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
	3938	#endif
	3939	if(ierr/=NF_NOERR) then
	3940	write(,) NF_STRERROR(ierr)
	3941	stop "getvarup"
	3942	endif
	3943	! write(,)'lecture dqtdy ok',dqtdy
	3944
	3945	#ifdef NC_DOUBLE
	3946	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
	3947	#else
	3948	ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
	3949	#endif
	3950	if(ierr/=NF_NOERR) then
	3951	write(,) NF_STRERROR(ierr)
	3952	stop "getvarup"
	3953	endif
	3954	! write(,)'lecture dqtdt ok',dqtdt
	3955
	3956	#ifdef NC_DOUBLE
	3957	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
	3958	#else
	3959	ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
	3960	#endif
	3961	if(ierr/=NF_NOERR) then
	3962	write(,) NF_STRERROR(ierr)
	3963	stop "getvarup"
	3964	endif
	3965	! write(,)'lecture thl_rad ok',thl_rad
	3966
	3967	return
	3968	end subroutine read_fire
[2096]	3969	!=====================================================================
[2126]	3970	subroutine read_dice(fich_dice,nlevel,ntime &
[2683]	3971	& ,zz,pres,t,qv,u,v,o3 &
[2126]	3972	& ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg &
	3973	& ,hadvt,hadvq,hadvu,hadvv,w,omega)
	3974
	3975	!program reading initial profils and forcings of the Dice case study
	3976
	3977
	3978	implicit none
	3979
	3980	#include "netcdf.inc"
[2683]	3981	#include "YOMCST.h"
[2126]	3982
	3983	integer ntime,nlevel
	3984	integer l,k
	3985	character*80 :: fich_dice
	3986	real*8 time(ntime)
	3987	real*8 zz(nlevel)
	3988
[2683]	3989	real*8 th(nlevel),pres(nlevel),t(nlevel)
[2126]	3990	real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
	3991	real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
	3992	real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
	3993	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
	3994	real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
[2683]	3995	real*8 pzero
[2126]	3996
	3997	integer nid, ierr
	3998	integer nbvar3d
	3999	parameter(nbvar3d=30)
	4000	integer var3didin(nbvar3d)
	4001
[2683]	4002	pzero=100000.
[2126]	4003	ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
	4004	if (ierr.NE.NF_NOERR) then
	4005	write(,) 'ERROR: Pb opening forcings nc file '
	4006	write(,) NF_STRERROR(ierr)
	4007	stop ""
	4008	endif
	4009
	4010
	4011	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
	4012	if(ierr/=NF_NOERR) then
	4013	write(,) NF_STRERROR(ierr)
	4014	stop 'height'
	4015	endif
	4016
	4017	ierr=NF_INQ_VARID(nid,"pf",var3didin(11))
	4018	if(ierr/=NF_NOERR) then
	4019	write(,) NF_STRERROR(ierr)
	4020	stop 'pf'
	4021	endif
	4022
	4023	ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
	4024	if(ierr/=NF_NOERR) then
	4025	write(,) NF_STRERROR(ierr)
	4026	stop 'theta'
	4027	endif
	4028
	4029	ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
	4030	if(ierr/=NF_NOERR) then
	4031	write(,) NF_STRERROR(ierr)
	4032	stop 'qv'
	4033	endif
	4034
	4035	ierr=NF_INQ_VARID(nid,"u",var3didin(14))
	4036	if(ierr/=NF_NOERR) then
	4037	write(,) NF_STRERROR(ierr)
	4038	stop 'u'
	4039	endif
	4040
	4041	ierr=NF_INQ_VARID(nid,"v",var3didin(15))
	4042	if(ierr/=NF_NOERR) then
	4043	write(,) NF_STRERROR(ierr)
	4044	stop 'v'
	4045	endif
	4046
	4047	ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
	4048	if(ierr/=NF_NOERR) then
	4049	write(,) NF_STRERROR(ierr)
	4050	stop 'o3'
	4051	endif
	4052
	4053	ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
	4054	if(ierr/=NF_NOERR) then
	4055	write(,) NF_STRERROR(ierr)
	4056	stop 'shf'
	4057	endif
	4058
	4059	ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
	4060	if(ierr/=NF_NOERR) then
	4061	write(,) NF_STRERROR(ierr)
	4062	stop 'lhf'
	4063	endif
	4064
	4065	ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
	4066	if(ierr/=NF_NOERR) then
	4067	write(,) NF_STRERROR(ierr)
	4068	stop 'lwup'
	4069	endif
	4070
	4071	ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
	4072	if(ierr/=NF_NOERR) then
	4073	write(,) NF_STRERROR(ierr)
	4074	stop 'dqtdx'
	4075	endif
	4076
	4077	ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
	4078	if(ierr/=NF_NOERR) then
	4079	write(,) NF_STRERROR(ierr)
	4080	stop 'Tg'
	4081	endif
	4082
	4083	ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
	4084	if(ierr/=NF_NOERR) then
	4085	write(,) NF_STRERROR(ierr)
	4086	stop 'ustar'
	4087	endif
	4088
	4089	ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
	4090	if(ierr/=NF_NOERR) then
	4091	write(,) NF_STRERROR(ierr)
	4092	stop 'psurf'
	4093	endif
	4094
	4095	ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
	4096	if(ierr/=NF_NOERR) then
	4097	write(,) NF_STRERROR(ierr)
	4098	stop 'Ug'
	4099	endif
	4100
	4101	ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
	4102	if(ierr/=NF_NOERR) then
	4103	write(,) NF_STRERROR(ierr)
	4104	stop 'Vg'
	4105	endif
	4106
	4107	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
	4108	if(ierr/=NF_NOERR) then
	4109	write(,) NF_STRERROR(ierr)
	4110	stop 'hadvT'
	4111	endif
	4112
	4113	ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
	4114	if(ierr/=NF_NOERR) then
	4115	write(,) NF_STRERROR(ierr)
	4116	stop 'hadvq'
	4117	endif
	4118
	4119	ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
	4120	if(ierr/=NF_NOERR) then
	4121	write(,) NF_STRERROR(ierr)
	4122	stop 'hadvu'
	4123	endif
	4124
	4125	ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
	4126	if(ierr/=NF_NOERR) then
	4127	write(,) NF_STRERROR(ierr)
	4128	stop 'hadvv'
	4129	endif
	4130
	4131	ierr=NF_INQ_VARID(nid,"w",var3didin(21))
	4132	if(ierr/=NF_NOERR) then
	4133	write(,) NF_STRERROR(ierr)
	4134	stop 'w'
	4135	endif
	4136
	4137	ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
	4138	if(ierr/=NF_NOERR) then
	4139	write(,) NF_STRERROR(ierr)
	4140	stop 'omega'
	4141	endif
	4142	!dimensions lecture
	4143	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
	4144
	4145	#ifdef NC_DOUBLE
	4146	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
	4147	#else
	4148	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
	4149	#endif
	4150	if(ierr/=NF_NOERR) then
	4151	write(,) NF_STRERROR(ierr)
	4152	stop "getvarup"
	4153	endif
	4154	! write(,)'lecture zz ok',zz
	4155
	4156	#ifdef NC_DOUBLE
	4157	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
	4158	#else
	4159	ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
	4160	#endif
	4161	if(ierr/=NF_NOERR) then
	4162	write(,) NF_STRERROR(ierr)
	4163	stop "getvarup"
	4164	endif
	4165	! write(,)'lecture pres ok',pres
	4166
	4167	#ifdef NC_DOUBLE
	4168	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
	4169	#else
	4170	ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
	4171	#endif
	4172	if(ierr/=NF_NOERR) then
	4173	write(,) NF_STRERROR(ierr)
	4174	stop "getvarup"
	4175	endif
	4176	! write(,)'lecture th ok',th
[2683]	4177	do k=1,nlevel
	4178	t(k)=th(k)(pres(k)/pzero)*rkappa
	4179	enddo
[2126]	4180
	4181	#ifdef NC_DOUBLE
	4182	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
	4183	#else
	4184	ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
	4185	#endif
	4186	if(ierr/=NF_NOERR) then
	4187	write(,) NF_STRERROR(ierr)
	4188	stop "getvarup"
	4189	endif
	4190	! write(,)'lecture qv ok',qv
	4191
	4192	#ifdef NC_DOUBLE
	4193	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
	4194	#else
	4195	ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
	4196	#endif
	4197	if(ierr/=NF_NOERR) then
	4198	write(,) NF_STRERROR(ierr)
	4199	stop "getvarup"
	4200	endif
	4201	! write(,)'lecture u ok',u
	4202
	4203	#ifdef NC_DOUBLE
	4204	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
	4205	#else
	4206	ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
	4207	#endif
	4208	if(ierr/=NF_NOERR) then
	4209	write(,) NF_STRERROR(ierr)
	4210	stop "getvarup"
	4211	endif
	4212	! write(,)'lecture v ok',v
	4213
	4214	#ifdef NC_DOUBLE
	4215	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
	4216	#else
	4217	ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
	4218	#endif
	4219	if(ierr/=NF_NOERR) then
	4220	write(,) NF_STRERROR(ierr)
	4221	stop "getvarup"
	4222	endif
	4223	! write(,)'lecture o3 ok',o3
	4224
	4225	#ifdef NC_DOUBLE
	4226	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
	4227	#else
	4228	ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
	4229	#endif
	4230	if(ierr/=NF_NOERR) then
	4231	write(,) NF_STRERROR(ierr)
	4232	stop "getvarup"
	4233	endif
	4234	! write(,)'lecture shf ok',shf
	4235
	4236	#ifdef NC_DOUBLE
	4237	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
	4238	#else
	4239	ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
	4240	#endif
	4241	if(ierr/=NF_NOERR) then
	4242	write(,) NF_STRERROR(ierr)
	4243	stop "getvarup"
	4244	endif
	4245	! write(,)'lecture lhf ok',lhf
	4246
	4247	#ifdef NC_DOUBLE
	4248	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
	4249	#else
	4250	ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
	4251	#endif
	4252	if(ierr/=NF_NOERR) then
	4253	write(,) NF_STRERROR(ierr)
	4254	stop "getvarup"
	4255	endif
	4256	! write(,)'lecture lwup ok',lwup
	4257
	4258	#ifdef NC_DOUBLE
	4259	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
	4260	#else
	4261	ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
	4262	#endif
	4263	if(ierr/=NF_NOERR) then
	4264	write(,) NF_STRERROR(ierr)
	4265	stop "getvarup"
	4266	endif
	4267	! write(,)'lecture swup ok',swup
	4268
	4269	#ifdef NC_DOUBLE
	4270	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
	4271	#else
	4272	ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
	4273	#endif
	4274	if(ierr/=NF_NOERR) then
	4275	write(,) NF_STRERROR(ierr)
	4276	stop "getvarup"
	4277	endif
	4278	! write(,)'lecture tg ok',tg
	4279
	4280	#ifdef NC_DOUBLE
	4281	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
	4282	#else
	4283	ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
	4284	#endif
	4285	if(ierr/=NF_NOERR) then
	4286	write(,) NF_STRERROR(ierr)
	4287	stop "getvarup"
	4288	endif
	4289	! write(,)'lecture ustar ok',ustar
	4290
	4291	#ifdef NC_DOUBLE
	4292	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
	4293	#else
	4294	ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
	4295	#endif
	4296	if(ierr/=NF_NOERR) then
	4297	write(,) NF_STRERROR(ierr)
	4298	stop "getvarup"
	4299	endif
	4300	! write(,)'lecture psurf ok',psurf
	4301
	4302	#ifdef NC_DOUBLE
	4303	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
	4304	#else
	4305	ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
	4306	#endif
	4307	if(ierr/=NF_NOERR) then
	4308	write(,) NF_STRERROR(ierr)
	4309	stop "getvarup"
	4310	endif
	4311	! write(,)'lecture ug ok',ug
	4312
	4313	#ifdef NC_DOUBLE
	4314	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
	4315	#else
	4316	ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
	4317	#endif
	4318	if(ierr/=NF_NOERR) then
	4319	write(,) NF_STRERROR(ierr)
	4320	stop "getvarup"
	4321	endif
	4322	! write(,)'lecture vg ok',vg
	4323
	4324	#ifdef NC_DOUBLE
	4325	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
	4326	#else
	4327	ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
	4328	#endif
	4329	if(ierr/=NF_NOERR) then
	4330	write(,) NF_STRERROR(ierr)
	4331	stop "getvarup"
	4332	endif
	4333	! write(,)'lecture hadvt ok',hadvt
	4334
	4335	#ifdef NC_DOUBLE
	4336	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
	4337	#else
	4338	ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
	4339	#endif
	4340	if(ierr/=NF_NOERR) then
	4341	write(,) NF_STRERROR(ierr)
	4342	stop "getvarup"
	4343	endif
	4344	! write(,)'lecture hadvq ok',hadvq
	4345
	4346	#ifdef NC_DOUBLE
	4347	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
	4348	#else
	4349	ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
	4350	#endif
	4351	if(ierr/=NF_NOERR) then
	4352	write(,) NF_STRERROR(ierr)
	4353	stop "getvarup"
	4354	endif
	4355	! write(,)'lecture hadvu ok',hadvu
	4356
	4357	#ifdef NC_DOUBLE
	4358	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
	4359	#else
	4360	ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
	4361	#endif
	4362	if(ierr/=NF_NOERR) then
	4363	write(,) NF_STRERROR(ierr)
	4364	stop "getvarup"
	4365	endif
	4366	! write(,)'lecture hadvv ok',hadvv
	4367
	4368	#ifdef NC_DOUBLE
	4369	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
	4370	#else
	4371	ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
	4372	#endif
	4373	if(ierr/=NF_NOERR) then
	4374	write(,) NF_STRERROR(ierr)
	4375	stop "getvarup"
	4376	endif
	4377	! write(,)'lecture w ok',w
	4378
	4379	#ifdef NC_DOUBLE
	4380	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
	4381	#else
	4382	ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
	4383	#endif
	4384	if(ierr/=NF_NOERR) then
	4385	write(,) NF_STRERROR(ierr)
	4386	stop "getvarup"
	4387	endif
	4388	! write(,)'lecture omega ok',omega
	4389
	4390	return
	4391	end subroutine read_dice
	4392	!=====================================================================
[2672]	4393	subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol &
	4394	& ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
	4395
	4396	!program reading initial profils and forcings of the Gabls4 case study
	4397
	4398
	4399	implicit none
	4400
	4401	#include "netcdf.inc"
	4402
	4403	integer ntime,nlevel,nsol
	4404	integer l,k
	4405	character*80 :: fich_gabls4
	4406	real*8 time(ntime)
	4407
	4408	! ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees
	4409	! dans un ordre inverse par rapport a la convention LMDZ
	4410	! ==> il faut tout inverser (MPL 20141024)
	4411	! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
	4412	real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
	4413	real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
	4414	real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
	4415
	4416	real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
	4417	real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
	4418	real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
	4419
	4420	real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
	4421	real*8 tg(ntime)
	4422	integer nid, ierr
	4423	integer nbvar3d
	4424	parameter(nbvar3d=30)
	4425	integer var3didin(nbvar3d)
	4426
	4427	ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
	4428	if (ierr.NE.NF_NOERR) then
	4429	write(,) 'ERROR: Pb opening forcings nc file '
	4430	write(,) NF_STRERROR(ierr)
	4431	stop ""
	4432	endif
	4433
	4434
	4435	ierr=NF_INQ_VARID(nid,"height",var3didin(1))
	4436	if(ierr/=NF_NOERR) then
	4437	write(,) NF_STRERROR(ierr)
	4438	stop 'height'
	4439	endif
	4440
	4441	ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
	4442	if(ierr/=NF_NOERR) then
	4443	write(,) NF_STRERROR(ierr)
	4444	stop 'depth_sn'
	4445	endif
	4446
	4447	ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
	4448	if(ierr/=NF_NOERR) then
	4449	write(,) NF_STRERROR(ierr)
	4450	stop 'Ug'
	4451	endif
	4452
	4453	ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
	4454	if(ierr/=NF_NOERR) then
	4455	write(,) NF_STRERROR(ierr)
	4456	stop 'Vg'
	4457	endif
	4458	ierr=NF_INQ_VARID(nid,"pf",var3didin(5))
	4459	if(ierr/=NF_NOERR) then
	4460	write(,) NF_STRERROR(ierr)
	4461	stop 'pf'
	4462	endif
	4463
	4464	ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
	4465	if(ierr/=NF_NOERR) then
	4466	write(,) NF_STRERROR(ierr)
	4467	stop 'theta'
	4468	endif
	4469
	4470	ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
	4471	if(ierr/=NF_NOERR) then
	4472	write(,) NF_STRERROR(ierr)
	4473	stop 'tempe'
	4474	endif
	4475
	4476	ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
	4477	if(ierr/=NF_NOERR) then
	4478	write(,) NF_STRERROR(ierr)
	4479	stop 'qv'
	4480	endif
	4481
	4482	ierr=NF_INQ_VARID(nid,"u",var3didin(9))
	4483	if(ierr/=NF_NOERR) then
	4484	write(,) NF_STRERROR(ierr)
	4485	stop 'u'
	4486	endif
	4487
	4488	ierr=NF_INQ_VARID(nid,"v",var3didin(10))
	4489	if(ierr/=NF_NOERR) then
	4490	write(,) NF_STRERROR(ierr)
	4491	stop 'v'
	4492	endif
	4493
	4494	ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
	4495	if(ierr/=NF_NOERR) then
	4496	write(,) NF_STRERROR(ierr)
	4497	stop 'hadvt'
	4498	endif
	4499
	4500	ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
	4501	if(ierr/=NF_NOERR) then
	4502	write(,) NF_STRERROR(ierr)
	4503	stop 'hadvq'
	4504	endif
	4505
	4506	ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
	4507	if(ierr/=NF_NOERR) then
	4508	write(,) NF_STRERROR(ierr)
	4509	stop 'tsnow'
	4510	endif
	4511
	4512	ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
	4513	if(ierr/=NF_NOERR) then
	4514	write(,) NF_STRERROR(ierr)
	4515	stop 'snow_density'
	4516	endif
	4517
	4518	ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
	4519	if(ierr/=NF_NOERR) then
	4520	write(,) NF_STRERROR(ierr)
	4521	stop 'Tg'
	4522	endif
	4523
	4524
	4525	!dimensions lecture
	4526	! call catchaxis(nid,ntime,nlevel,time,z,ierr)
	4527
	4528	#ifdef NC_DOUBLE
	4529	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)
	4530	#else
	4531	ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)
	4532	#endif
	4533	if(ierr/=NF_NOERR) then
	4534	write(,) NF_STRERROR(ierr)
	4535	stop "getvarup"
	4536	endif
	4537
	4538	#ifdef NC_DOUBLE
	4539	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)
	4540	#else
	4541	ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)
	4542	#endif
	4543	if(ierr/=NF_NOERR) then
	4544	write(,) NF_STRERROR(ierr)
	4545	stop "getvarup"
	4546	endif
	4547
	4548	#ifdef NC_DOUBLE
	4549	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)
	4550	#else
	4551	ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)
	4552	#endif
	4553	if(ierr/=NF_NOERR) then
	4554	write(,) NF_STRERROR(ierr)
	4555	stop "getvarup"
	4556	endif
	4557
	4558	#ifdef NC_DOUBLE
	4559	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)
	4560	#else
	4561	ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)
	4562	#endif
	4563	if(ierr/=NF_NOERR) then
	4564	write(,) NF_STRERROR(ierr)
	4565	stop "getvarup"
	4566	endif
	4567
	4568	#ifdef NC_DOUBLE
	4569	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)
	4570	#else
	4571	ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)
	4572	#endif
	4573	if(ierr/=NF_NOERR) then
	4574	write(,) NF_STRERROR(ierr)
	4575	stop "getvarup"
	4576	endif
	4577
	4578	#ifdef NC_DOUBLE
	4579	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)
	4580	#else
	4581	ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)
	4582	#endif
	4583	if(ierr/=NF_NOERR) then
	4584	write(,) NF_STRERROR(ierr)
	4585	stop "getvarup"
	4586	endif
	4587
	4588	#ifdef NC_DOUBLE
	4589	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)
	4590	#else
	4591	ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)
	4592	#endif
	4593	if(ierr/=NF_NOERR) then
	4594	write(,) NF_STRERROR(ierr)
	4595	stop "getvarup"
	4596	endif
	4597
	4598	#ifdef NC_DOUBLE
	4599	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)
	4600	#else
	4601	ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)
	4602	#endif
	4603	if(ierr/=NF_NOERR) then
	4604	write(,) NF_STRERROR(ierr)
	4605	stop "getvarup"
	4606	endif
	4607
	4608	#ifdef NC_DOUBLE
	4609	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)
	4610	#else
	4611	ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)
	4612	#endif
	4613	if(ierr/=NF_NOERR) then
	4614	write(,) NF_STRERROR(ierr)
	4615	stop "getvarup"
	4616	endif
	4617
	4618	#ifdef NC_DOUBLE
	4619	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)
	4620	#else
	4621	ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)
	4622	#endif
	4623	if(ierr/=NF_NOERR) then
	4624	write(,) NF_STRERROR(ierr)
	4625	stop "getvarup"
	4626	endif
	4627
	4628	#ifdef NC_DOUBLE
	4629	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)
	4630	#else
	4631	ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)
	4632	#endif
	4633	if(ierr/=NF_NOERR) then
	4634	write(,) NF_STRERROR(ierr)
	4635	stop "getvarup"
	4636	endif
	4637
	4638	#ifdef NC_DOUBLE
	4639	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)
	4640	#else
	4641	ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)
	4642	#endif
	4643	if(ierr/=NF_NOERR) then
	4644	write(,) NF_STRERROR(ierr)
	4645	stop "getvarup"
	4646	endif
	4647
	4648	#ifdef NC_DOUBLE
	4649	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)
	4650	#else
	4651	ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)
	4652	#endif
	4653	if(ierr/=NF_NOERR) then
	4654	write(,) NF_STRERROR(ierr)
	4655	stop "getvarup"
	4656	endif
	4657
	4658	#ifdef NC_DOUBLE
	4659	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)
	4660	#else
	4661	ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)
	4662	#endif
	4663	if(ierr/=NF_NOERR) then
	4664	write(,) NF_STRERROR(ierr)
	4665	stop "getvarup"
	4666	endif
	4667
	4668	#ifdef NC_DOUBLE
	4669	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)
	4670	#else
	4671	ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)
	4672	#endif
	4673	if(ierr/=NF_NOERR) then
	4674	write(,) NF_STRERROR(ierr)
	4675	stop "getvarup"
	4676	endif
	4677
	4678	! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
	4679	do k=1,nlevel
	4680	zz(k)=zz_i(nlevel+1-k)
	4681	ug(k,:)=ug_i(nlevel+1-k,:)
	4682	vg(k,:)=vg_i(nlevel+1-k,:)
	4683	pf(k)=pf_i(nlevel+1-k)
	4684	print *,'pf=',pf(k)
	4685	th(k)=th_i(nlevel+1-k)
	4686	t(k)=t_i(nlevel+1-k)
	4687	qv(k)=qv_i(nlevel+1-k)
	4688	u(k)=u_i(nlevel+1-k)
	4689	v(k)=v_i(nlevel+1-k)
	4690	hadvt(k,:)=hadvt_i(nlevel+1-k,:)
	4691	hadvq(k,:)=hadvq_i(nlevel+1-k,:)
	4692	enddo
	4693	return
	4694	end subroutine read_gabls4
	4695	!=====================================================================
	4696
[2096]	4697	! Reads CIRC input files
	4698
	4699	SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
	4700
	4701	parameter (ncm_1=49180)
	4702	#include "YOMCST.h"
	4703
	4704	real albsfc(ncm_1), albsfc_w(ncm_1)
	4705	real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
	4706	reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
	4707	real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
	4708	real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
	4709	real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
	4710	real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
	4711	o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
	4712	! za= zenital angle
	4713	! sza= cosinus angle zenital
	4714	real wavn(ncm_1), ssf(ncm_1),za,sza
	4715	integer nlev
	4716
	4717
	4718	! Open the files
	4719
	4720	open (11, file='Tsfc_sza_nlev_case.txt', status='old')
	4721	open (12, file='level_input_case.txt', status='old')
	4722	open (13, file='layer_input_case.txt', status='old')
	4723	open (14, file='aerosol_input_case.txt', status='old')
	4724	open (15, file='cloud_input_case.txt', status='old')
	4725	open (16, file='sfcalbedo_input_case.txt', status='old')
	4726
	4727	! Read scalar information
	4728	do iskip=1,5
	4729	read (11, *)
	4730	enddo
	4731	read (11, '(i8)') nlev
	4732	read (11, '(f10.2)') tsfc
	4733	read (11, '(f10.2)') za
	4734	read (11, '(f10.4)') sw_dn_toa
	4735	sza=cos(za/180.*RPI)
	4736	print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
	4737	close(11)
	4738
	4739	! Read level information
	4740	read (12, *)
	4741	do il=1,nlev
	4742	read (12, 302) ilev, z(il), p(il), t(il)
	4743	z(il)=z(il)*1000. ! z donne en km
	4744	p(il)=p(il)*100. ! p donne en mb
	4745	enddo
	4746	302 format (i8, f8.3, 2f9.2)
	4747	close(12)
	4748
	4749	! Read layer information (midpoint values)
	4750	do iskip=1,3
	4751	read (13, *)
	4752	enddo
	4753	do il=1,nlev-1
	4754	read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
	4755	n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
	4756	f11(il),f12(il)
	4757	pm(il)=pm(il)*100.
	4758	enddo
	4759	303 format (i8, 2f9.2, 10(2x,e13.7))
	4760	close(13)
	4761
	4762	! Read aerosol layer information
	4763	do iskip=1,3
	4764	read (14, *)
	4765	enddo
	4766	read (14, '(f10.2)') aer_alpha
	4767	read (14, *)
	4768	read (14, *)
	4769	do il=1,nlev-1
	4770	read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
	4771	enddo
	4772	304 format (i8, f9.5, 2f8.3)
	4773	close(14)
	4774
	4775	! Read cloud information
	4776	do iskip=1,3
	4777	read (15, *)
	4778	enddo
	4779	do il=1,nlev-1
	4780	read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
	4781	lwp(il)=lwp(il)/1000. ! lwp donne en g/kg
	4782	iwp(il)=iwp(il)/1000. ! iwp donne en g/kg
	4783	reliq(il)=reliq(il)/1000000. ! reliq donne en microns
	4784	reice(il)=reice(il)/1000000. ! reice donne en microns
	4785	enddo
	4786	305 format (i8, f8.3, 4f9.2)
	4787	close(15)
	4788
	4789	! Read surface albedo (weighted & unweighted) and spectral solar irradiance
	4790	do iskip=1,6
	4791	read (16, *)
	4792	enddo
	4793	do icm_1=1,ncm_1
	4794	read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
	4795	enddo
	4796	306 format(f10.1, 2f12.5, f14.8)
	4797	close(16)
	4798
	4799	return
	4800	end subroutine read_circ
	4801	!=====================================================================
	4802	! Reads RTMIP input files
	4803
	4804	SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
	4805
	4806	#include "YOMCST.h"
	4807
	4808	real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
	4809	real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
	4810	integer nlev
	4811
	4812
	4813	! Open the files
	4814
	4815	open (11, file='low_resolution_profile.txt', status='old')
	4816
	4817	! Read level information
	4818	read (11, *)
	4819	do il=1,nlev_rtmip
	4820	read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
	4821	enddo
	4822	do il=1,nlev_rtmip
	4823	play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb
	4824	temp(il)=t(nlev_rtmip-il+1)
	4825	ovap(il)=h2o(nlev_rtmip-il+1)
	4826	oz(il)=o3(nlev_rtmip-il+1)
	4827	enddo
	4828	do il=1,39
	4829	plev(il)=play(il)+(play(il+1)-play(il))/2.
	4830	print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
	4831	enddo
	4832	plev(41)=101300.
	4833	302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
	4834	close(12)
	4835
	4836	return
	4837	end subroutine read_rtmip
	4838	!=====================================================================
	4839
[2181]	4840	! Subroutines for nudging
	4841
	4842	Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
	4843	! ========================================================
	4844	USE dimphy
	4845
	4846	implicit none
	4847
	4848	! ========================================================
	4849	REAL paprs(klon,klevp1)
	4850	REAL pplay(klon,klev)
	4851	!
	4852	! Variables d'etat
	4853	REAL t(klon,klev)
	4854	REAL q(klon,klev)
	4855	!
	4856	! Profiles cible
	4857	REAL t_targ(klon,klev)
	4858	REAL rh_targ(klon,klev)
	4859	!
	4860	INTEGER k,i
	4861	REAL zx_qs
	4862
	4863	! Declaration des constantes et des fonctions thermodynamiques
	4864	!
	4865	include "YOMCST.h"
	4866	include "YOETHF.h"
	4867	!
	4868	! ----------------------------------------
	4869	! Statement functions
	4870	include "FCTTRE.h"
	4871	! ----------------------------------------
	4872	!
	4873	DO k = 1,klev
	4874	DO i = 1,klon
	4875	t_targ(i,k) = t(i,k)
	4876	IF (t(i,k).LT.RTT) THEN
	4877	zx_qs = qsats(t(i,k))/(pplay(i,k))
	4878	ELSE
	4879	zx_qs = qsatl(t(i,k))/(pplay(i,k))
	4880	ENDIF
	4881	rh_targ(i,k) = q(i,k)/zx_qs
	4882	ENDDO
	4883	ENDDO
	4884	print *, 't_targ',t_targ
	4885	print *, 'rh_targ',rh_targ
	4886	!
	4887	!
	4888	RETURN
	4889	END
	4890
	4891	Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
	4892	! ========================================================
	4893	USE dimphy
	4894
	4895	implicit none
	4896
	4897	! ========================================================
	4898	REAL paprs(klon,klevp1)
	4899	REAL pplay(klon,klev)
	4900	!
	4901	! Variables d'etat
	4902	REAL u(klon,klev)
	4903	REAL v(klon,klev)
	4904	!
	4905	! Profiles cible
	4906	REAL u_targ(klon,klev)
	4907	REAL v_targ(klon,klev)
	4908	!
	4909	INTEGER k,i
	4910	!
	4911	DO k = 1,klev
	4912	DO i = 1,klon
	4913	u_targ(i,k) = u(i,k)
	4914	v_targ(i,k) = v(i,k)
	4915	ENDDO
	4916	ENDDO
	4917	print *, 'u_targ',u_targ
	4918	print *, 'v_targ',v_targ
	4919	!
	4920	!
	4921	RETURN
	4922	END
	4923
	4924	Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q, &
	4925	& d_t,d_q)
	4926	! ========================================================
	4927	USE dimphy
	4928
	4929	implicit none
	4930
	4931	! ========================================================
	4932	REAL dtime
	4933	REAL paprs(klon,klevp1)
	4934	REAL pplay(klon,klev)
	4935	!
	4936	! Variables d'etat
	4937	REAL t(klon,klev)
	4938	REAL q(klon,klev)
	4939	!
	4940	! Tendances
	4941	REAL d_t(klon,klev)
	4942	REAL d_q(klon,klev)
	4943	!
	4944	! Profiles cible
	4945	REAL t_targ(klon,klev)
	4946	REAL rh_targ(klon,klev)
	4947	!
	4948	! Temps de relaxation
	4949	REAL tau
	4950	!c DATA tau /3600./
	4951	!! DATA tau /5400./
	4952	DATA tau /1800./
	4953	!
	4954	INTEGER k,i
[2455]	4955	REAL zx_qs, rh, tnew, d_rh, rhnew
[2181]	4956
	4957	! Declaration des constantes et des fonctions thermodynamiques
	4958	!
	4959	include "YOMCST.h"
	4960	include "YOETHF.h"
	4961	!
	4962	! ----------------------------------------
	4963	! Statement functions
	4964	include "FCTTRE.h"
	4965	! ----------------------------------------
	4966	!
	4967	print *,'dtime, tau ',dtime,tau
	4968	print *, 't_targ',t_targ
	4969	print *, 'rh_targ',rh_targ
	4970	print *,'temp ',t
	4971	print *,'hum ',q
[2455]	4972	!
[2181]	4973	DO k = 1,klev
	4974	DO i = 1,klon
[2455]	4975	IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
[2181]	4976	IF (t(i,k).LT.RTT) THEN
	4977	zx_qs = qsats(t(i,k))/(pplay(i,k))
	4978	ELSE
	4979	zx_qs = qsatl(t(i,k))/(pplay(i,k))
	4980	ENDIF
	4981	rh = q(i,k)/zx_qs
	4982	!
	4983	d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
	4984	d_rh = 1./tau*(rh_targ(i,k)-rh)
	4985	!
[2455]	4986	tnew = t(i,k)+d_t(i,k)*dtime
	4987	!jyg<
	4988	! Formule pour q :
	4989	! d_q = (1/tau) [rh_targ*qsat(T_new) - q]
	4990	!
	4991	! Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
[2716]	4992	! qui n'etait pas correcte.
[2455]	4993	!
[2181]	4994	IF (tnew.LT.RTT) THEN
	4995	zx_qs = qsats(tnew)/(pplay(i,k))
	4996	ELSE
	4997	zx_qs = qsatl(tnew)/(pplay(i,k))
	4998	ENDIF
[2455]	4999	!! d_q(i,k) = d_q(i,k) + d_rh*zx_qs
	5000	d_q(i,k) = d_q(i,k) + (1./tau)(rh_targ(i,k)zx_qs - q(i,k))
	5001	rhnew = (q(i,k)+d_q(i,k)*dtime)/zx_qs
[2181]	5002	!
[2455]	5003	print *,' k,d_t,rh,d_rh,rhnew,d_q ', &
	5004	k,d_t(i,k),rh,d_rh,rhnew,d_q(i,k)
	5005	ENDIF
[2181]	5006	!
	5007	ENDDO
	5008	ENDDO
	5009	!
	5010	RETURN
	5011	END
	5012
	5013	Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v, &
	5014	& d_u,d_v)
	5015	! ========================================================
	5016	USE dimphy
	5017
	5018	implicit none
	5019
	5020	! ========================================================
	5021	REAL dtime
	5022	REAL paprs(klon,klevp1)
	5023	REAL pplay(klon,klev)
	5024	!
	5025	! Variables d'etat
	5026	REAL u(klon,klev)
	5027	REAL v(klon,klev)
	5028	!
	5029	! Tendances
	5030	REAL d_u(klon,klev)
	5031	REAL d_v(klon,klev)
	5032	!
	5033	! Profiles cible
	5034	REAL u_targ(klon,klev)
	5035	REAL v_targ(klon,klev)
	5036	!
	5037	! Temps de relaxation
	5038	REAL tau
	5039	!c DATA tau /3600./
[2933]	5040	! DATA tau /5400./
	5041	DATA tau /43200./
[2181]	5042	!
	5043	INTEGER k,i
	5044
	5045	!
	5046	print *,'dtime, tau ',dtime,tau
	5047	print *, 'u_targ',u_targ
	5048	print *, 'v_targ',v_targ
	5049	print *,'zonal velocity ',u
	5050	print *,'meridional velocity ',v
	5051	DO k = 1,klev
	5052	DO i = 1,klon
[2933]	5053	!CR: nudging everywhere
	5054	! IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
[2181]	5055	!
	5056	d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
	5057	d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
	5058	!
	5059	print *,' k,u,d_u,v,d_v ', &
	5060	k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
[2933]	5061	! ENDIF
[2181]	5062	!
	5063	ENDDO
	5064	ENDDO
	5065	!
	5066	RETURN
	5067	END
	5068
[2716]	5069	!=====================================================================
	5070	SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas &
	5071	& ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas &
	5072	& ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas &
	5073	& ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas &
	5074	& ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &
	5075	& ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
	5076	& ,dth_prof_cas,hth_prof_cas,vth_prof_cas &
	5077	!
	5078	& ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas &
	5079	& ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas &
	5080	& ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas &
	5081	& ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &
	5082	& ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas &
	5083	& ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
	5084
	5085	implicit none
	5086
[2990]	5087	#include "YOMCST.h"
[2716]	5088	#include "dimensions.h"
[2683]	5089
[2716]	5090	!-------------------------------------------------------------------------
	5091	! Vertical interpolation of generic case forcing data onto mod_casel levels
	5092	!-------------------------------------------------------------------------
	5093
	5094	integer nlevmax
	5095	parameter (nlevmax=41)
	5096	integer nlev_cas,mxcalc
	5097	! real play(llm), plev_prof(nlevmax)
	5098	! real t_prof(nlevmax),q_prof(nlevmax)
	5099	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	5100	! real ht_prof(nlevmax),vt_prof(nlevmax)
	5101	! real hq_prof(nlevmax),vq_prof(nlevmax)
	5102
	5103	real play(llm), plev_prof_cas(nlev_cas)
	5104	real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
	5105	real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
	5106	real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
	5107	real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
	5108	real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
	5109	real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
	5110	real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
	5111	real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
	5112	real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
	5113
	5114	real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
	5115	real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
	5116	real u_mod_cas(llm),v_mod_cas(llm)
	5117	real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
	5118	real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
	5119	real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
	5120	real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
	5121	real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
	5122	real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
	5123
	5124	integer l,k,k1,k2
	5125	real frac,frac1,frac2,fact
	5126
	5127	do l = 1, llm
	5128	print *,'debut interp2, play=',l,play(l)
	5129	enddo
	5130	! do l = 1, nlev_cas
	5131	! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
	5132	! enddo
	5133
	5134	do l = 1, llm
	5135
	5136	if (play(l).ge.plev_prof_cas(nlev_cas)) then
	5137
	5138	mxcalc=l
	5139	print *,'debut interp2, mxcalc=',mxcalc
	5140	k1=0
	5141	k2=0
	5142
	5143	if (play(l).le.plev_prof_cas(1)) then
	5144
	5145	do k = 1, nlev_cas-1
	5146	if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
	5147	k1=k
	5148	k2=k+1
	5149	endif
	5150	enddo
	5151
	5152	if (k1.eq.0 .or. k2.eq.0) then
	5153	write(,) 'PB! k1, k2 = ',k1,k2
	5154	write(,) 'l,play(l) = ',l,play(l)/100
	5155	do k = 1, nlev_cas-1
	5156	write(,) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
	5157	enddo
	5158	endif
	5159
	5160	frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
	5161	t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
	5162	theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
[3008]	5163	if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)(play(l)/100000.)*(RD/RCPD)
[2716]	5164	thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
	5165	thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
	5166	qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
	5167	ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
	5168	qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
	5169	u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
	5170	v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
	5171	ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
	5172	vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
	5173	w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
	5174	omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
	5175	du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
	5176	hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
	5177	vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
	5178	dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
	5179	hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
	5180	vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
	5181	dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
	5182	ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
	5183	vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
	5184	dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
	5185	hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
	5186	vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
	5187	dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
	5188	hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
	5189	vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
[2933]	5190	dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
[2716]	5191
	5192	else !play>plev_prof_cas(1)
	5193
	5194	k1=1
	5195	k2=2
	5196	print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
	5197	frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
	5198	frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
	5199	t_mod_cas(l)= frac1t_prof_cas(k1) - frac2t_prof_cas(k2)
	5200	theta_mod_cas(l)= frac1th_prof_cas(k1) - frac2th_prof_cas(k2)
[3008]	5201	if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)(play(l)/100000.)*(RD/RCPD)
[2716]	5202	thv_mod_cas(l)= frac1thv_prof_cas(k1) - frac2thv_prof_cas(k2)
	5203	thl_mod_cas(l)= frac1thl_prof_cas(k1) - frac2thl_prof_cas(k2)
	5204	qv_mod_cas(l)= frac1qv_prof_cas(k1) - frac2qv_prof_cas(k2)
	5205	ql_mod_cas(l)= frac1ql_prof_cas(k1) - frac2ql_prof_cas(k2)
	5206	qi_mod_cas(l)= frac1qi_prof_cas(k1) - frac2qi_prof_cas(k2)
	5207	u_mod_cas(l)= frac1u_prof_cas(k1) - frac2u_prof_cas(k2)
	5208	v_mod_cas(l)= frac1v_prof_cas(k1) - frac2v_prof_cas(k2)
	5209	ug_mod_cas(l)= frac1ug_prof_cas(k1) - frac2ug_prof_cas(k2)
	5210	vg_mod_cas(l)= frac1vg_prof_cas(k1) - frac2vg_prof_cas(k2)
	5211	w_mod_cas(l)= frac1vitw_prof_cas(k1) - frac2vitw_prof_cas(k2)
	5212	omega_mod_cas(l)= frac1omega_prof_cas(k1) - frac2omega_prof_cas(k2)
	5213	du_mod_cas(l)= frac1du_prof_cas(k1) - frac2du_prof_cas(k2)
	5214	hu_mod_cas(l)= frac1hu_prof_cas(k1) - frac2hu_prof_cas(k2)
	5215	vu_mod_cas(l)= frac1vu_prof_cas(k1) - frac2vu_prof_cas(k2)
	5216	dv_mod_cas(l)= frac1dv_prof_cas(k1) - frac2dv_prof_cas(k2)
	5217	hv_mod_cas(l)= frac1hv_prof_cas(k1) - frac2hv_prof_cas(k2)
	5218	vv_mod_cas(l)= frac1vv_prof_cas(k1) - frac2vv_prof_cas(k2)
	5219	dt_mod_cas(l)= frac1dt_prof_cas(k1) - frac2dt_prof_cas(k2)
	5220	ht_mod_cas(l)= frac1ht_prof_cas(k1) - frac2ht_prof_cas(k2)
	5221	vt_mod_cas(l)= frac1vt_prof_cas(k1) - frac2vt_prof_cas(k2)
	5222	dth_mod_cas(l)= frac1dth_prof_cas(k1) - frac2dth_prof_cas(k2)
	5223	hth_mod_cas(l)= frac1hth_prof_cas(k1) - frac2hth_prof_cas(k2)
	5224	vth_mod_cas(l)= frac1vth_prof_cas(k1) - frac2vth_prof_cas(k2)
	5225	dq_mod_cas(l)= frac1dq_prof_cas(k1) - frac2dq_prof_cas(k2)
	5226	hq_mod_cas(l)= frac1hq_prof_cas(k1) - frac2hq_prof_cas(k2)
	5227	vq_mod_cas(l)= frac1vq_prof_cas(k1) - frac2vq_prof_cas(k2)
[2933]	5228	dtrad_mod_cas(l)= frac1dtrad_prof_cas(k1) - frac2dtrad_prof_cas(k2)
[2716]	5229
	5230	endif ! play.le.plev_prof_cas(1)
	5231
	5232	else ! above max altitude of forcing file
	5233
	5234	!jyg
	5235	fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
	5236	fact = max(fact,0.) !jyg
	5237	fact = exp(-fact) !jyg
	5238	t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg
	5239	theta_mod_cas(l)= th_prof_cas(nlev_cas) !jyg
	5240	thv_mod_cas(l)= thv_prof_cas(nlev_cas) !jyg
	5241	thl_mod_cas(l)= thl_prof_cas(nlev_cas) !jyg
	5242	qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact !jyg
	5243	ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact !jyg
	5244	qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact !jyg
	5245	u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg
	5246	v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg
	5247	ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg
	5248	vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg
	5249	w_mod_cas(l)= 0.0 !jyg
[2990]	5250	omega_mod_cas(l)= 0.0 !jyg
[2716]	5251	du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
	5252	hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg
	5253	vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg
	5254	dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
	5255	hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg
	5256	vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg
	5257	dt_mod_cas(l)= dt_prof_cas(nlev_cas)
	5258	ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg
	5259	vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg
	5260	dth_mod_cas(l)= dth_prof_cas(nlev_cas)
	5261	hth_mod_cas(l)= hth_prof_cas(nlev_cas) !jyg
	5262	vth_mod_cas(l)= vth_prof_cas(nlev_cas) !jyg
	5263	dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
	5264	hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg
	5265	vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg
[2933]	5266	dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg
[2716]	5267
	5268	endif ! play
	5269
	5270	enddo ! l
	5271
	5272	return
	5273	end
	5274	!*****************************************************************************
	5275
	5276
[3008]	5277
	5278

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

Download in other formats: