Context Navigation

source: LMDZ5/trunk/libf/phylmd/1DUTILS.h @ 2103

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

bug dans la lecture du profil initial de Qv

MPL, PM

File size: 115.9 KB

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

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

Download in other formats: