Context Navigation

source: LMDZ5/trunk/libf/phy1d/1DUTILS.h @ 2007

Last change on this file since 2007 was 1979, checked in by fhourdin, 11 years ago
Suprresion de inigrads de 1DUTILS.h
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 108.1 KB

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

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

Download in other formats: