Context Navigation

1DUTILS.h_with_writelim @ 1618

Last change on this file since 1618 was 1608, checked in by lguez, 13 years ago

In "phy1d/*", replaced obsolete calls to "lnblnk" by calls to "trim".

In "lmdz1d.F", use collector module "ioipsl" rather than specific
module "calendar", which is an internal module of IOIPSL. Removed line
containing only "#" (causes compilation error). Bug fix in call to
"init_phys_lmdz".

File size: 80.5 KB

Rev	Line
[1607]	1	!
	2	! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
	3	!
	4	c
	5	c
	6	SUBROUTINE conf_unicol( tapedef )
	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 Arguments :
	19	c
	20	c tapedef :
	21
	22	INTEGER tapedef
	23	c
	24	c Declarations :
	25	c --------------
	26
	27	#include "compar1d.h"
	28	#include "flux_arp.h"
	29	#include "fcg_gcssold.h"
	30	#include "fcg_racmo.h"
	31	#include "iniprint.h"
	32	c
	33	c
	34	c local:
	35	c ------
	36
	37	c CHARACTER ch172,ch272,ch372,ch412
	38
	39	c
	40	c -------------------------------------------------------------------
	41	c
	42	c ......... Initilisation parametres du lmdz1D ..........
	43	c
	44	c---------------------------------------------------------------------
	45	c initialisations:
	46	c ----------------
	47
	48	!Config Key = lunout
	49	!Config Desc = unite de fichier pour les impressions
	50	!Config Def = 6
	51	!Config Help = unite de fichier pour les impressions
	52	!Config (defaut sortie standard = 6)
	53	lunout=6
	54	CALL getin('lunout', lunout)
	55	IF (lunout /= 5 .and. lunout /= 6) THEN
	56	OPEN(lunout,FILE='lmdz.out')
	57	ENDIF
	58
	59	!Config Key = prt_level
	60	!Config Desc = niveau d'impressions de débogage
	61	!Config Def = 0
	62	!Config Help = Niveau d'impression pour le débogage
	63	!Config (0 = minimum d'impression)
	64	c prt_level = 0
	65	c CALL getin('prt_level',prt_level)
	66
	67	c-----------------------------------------------------------------------
	68	c Parametres de controle du run:
	69	c-----------------------------------------------------------------------
	70
	71	!Config Key = restart
	72	!Config Desc = on repart des startphy et start1dyn
	73	!Config Def = false
	74	!Config Help = les fichiers restart doivent etre renomme en start
	75	restart = .FALSE.
	76	CALL getin('restart',restart)
	77
	78	!Config Key = forcing_type
	79	!Config Desc = defines the way the SCM is forced:
	80	!Config Def = 0
	81	!!Config Help = 0 ==> forcing_les = .true.
	82	! initial profiles from file prof.inp.001
	83	! no forcing by LS convergence ;
	84	! surface temperature imposed ;
	85	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
	86	! = 1 ==> forcing_radconv = .true.
	87	! idem forcing_type = 0, but the imposed radiative cooling
	88	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
	89	! then there is no radiative cooling at all)
	90	! = 2 ==> forcing_toga = .true.
	91	! initial profiles from TOGA-COARE IFA files
	92	! LS convergence and SST imposed from TOGA-COARE IFA files
	93	! = 3 ==> forcing_GCM2SCM = .true.
	94	! initial profiles from the GCM output
	95	! LS convergence imposed from the GCM output
	96	! = 4 ==> forcing_twpi = .true.
	97	! initial profiles from TWPICE nc files
	98	! LS convergence and SST imposed from TWPICE nc files
	99	! = 5 ==> forcing_rico = .true.
	100	! initial profiles from RICO idealized
	101	! LS convergence imposed from RICO (cst)
	102	! = 40 ==> forcing_GCSSold = .true.
	103	! initial profile from GCSS file
	104	! LS convergence imposed from GCSS file
	105	!
	106	forcing_type = 0
	107	CALL getin('forcing_type',forcing_type)
	108	imp_fcg_gcssold = .false.
	109	ts_fcg_gcssold = .false.
	110	Tp_fcg_gcssold = .false.
	111	Tp_ini_gcssold = .false.
	112	xTurb_fcg_gcssold = .false.
	113	IF (forcing_type .eq.40) THEN
	114	CALL getin('imp_fcg',imp_fcg_gcssold)
	115	CALL getin('ts_fcg',ts_fcg_gcssold)
	116	CALL getin('tp_fcg',Tp_fcg_gcssold)
	117	CALL getin('tp_ini',Tp_ini_gcssold)
	118	CALL getin('turb_fcg',xTurb_fcg_gcssold)
	119	ENDIF
	120
	121	!Config Key = ok_flux_surf
	122	!Config Desc = forcage ou non par les flux de surface
	123	!Config Def = false
	124	!Config Help = forcage ou non par les flux de surface
	125	ok_flux_surf = .FALSE.
	126	CALL getin('ok_flux_surf',ok_flux_surf)
	127
	128	!Config Key = time_ini
	129	!Config Desc = meaningless in this case
	130	!Config Def = 0.
	131	!Config Help =
	132	tsurf = 0.
	133	CALL getin('time_ini',time_ini)
	134
	135	!Config Key = rlat et rlon
	136	!Config Desc = latitude et longitude
	137	!Config Def = 0.0 0.0
	138	!Config Help = fixe la position de la colonne
	139	xlat = 0.
	140	xlon = 0.
	141	CALL getin('rlat',xlat)
	142	CALL getin('rlon',xlon)
	143
	144	!Config Key = airephy
	145	!Config Desc = Grid cell area
	146	!Config Def = 1.e11
	147	!Config Help =
	148	airefi = 1.e11
	149	CALL getin('airephy',airefi)
	150
	151	!Config Key = nat_surf
	152	!Config Desc = surface type
	153	!Config Def = 0 (ocean)
	154	!Config Help = 0=ocean,1=land,2=glacier,3=banquise
	155	nat_surf = 0.
	156	CALL getin('nat_surf',nat_surf)
	157
	158	!Config Key = tsurf
	159	!Config Desc = surface temperature
	160	!Config Def = 290.
	161	!Config Help = not used if type_ts_forcing=1 in lmdz1d.F
	162	tsurf = 290.
	163	CALL getin('tsurf',tsurf)
	164
	165	!Config Key = psurf
	166	!Config Desc = surface pressure
	167	!Config Def = 102400.
	168	!Config Help =
	169	psurf = 102400.
	170	CALL getin('psurf',psurf)
	171
	172	!Config Key = zsurf
	173	!Config Desc = surface altitude
	174	!Config Def = 0.
	175	!Config Help =
	176	zsurf = 0.
	177	CALL getin('zsurf',zsurf)
	178
	179	!Config Key = rugos
	180	!Config Desc = coefficient de frottement
	181	!Config Def = 0.0001
	182	!Config Help = calcul du Cdrag
	183	rugos = 0.0001
	184	CALL getin('rugos',rugos)
	185
	186	!Config Key = wtsurf et wqsurf
	187	!Config Desc = ???
	188	!Config Def = 0.0 0.0
	189	!Config Help =
	190	wtsurf = 0.0
	191	wqsurf = 0.0
	192	CALL getin('wtsurf',wtsurf)
	193	CALL getin('wqsurf',wqsurf)
	194
	195	!Config Key = albedo
	196	!Config Desc = albedo
	197	!Config Def = 0.09
	198	!Config Help =
	199	albedo = 0.09
	200	CALL getin('albedo',albedo)
	201
	202	!Config Key = agesno
	203	!Config Desc = age de la neige
	204	!Config Def = 30.0
	205	!Config Help =
	206	xagesno = 30.0
	207	CALL getin('agesno',xagesno)
	208
	209	!Config Key = restart_runoff
	210	!Config Desc = age de la neige
	211	!Config Def = 30.0
	212	!Config Help =
	213	restart_runoff = 0.0
	214	CALL getin('restart_runoff',restart_runoff)
	215
	216	!Config Key = qsolinp
	217	!Config Desc = initial bucket water content (kg/m2) when land (5std)
	218	!Config Def = 30.0
	219	!Config Help =
	220	qsolinp = 1.
	221	CALL getin('qsolinp',qsolinp)
	222
	223	!Config Key = zpicinp
	224	!Config Desc = denivellation orographie
	225	!Config Def = 300.
	226	!Config Help = input brise
	227	zpicinp = 300.
	228	CALL getin('zpicinp',zpicinp)
	229
	230
	231	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	232	write(lunout,*)' Configuration des parametres du gcm1D: '
	233	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	234	write(lunout,*)' restart = ', restart
	235	write(lunout,*)' forcing_type = ', forcing_type
	236	write(lunout,*)' time_ini = ', time_ini
	237	write(lunout,*)' rlat = ', xlat
	238	write(lunout,*)' rlon = ', xlon
	239	write(lunout,*)' airephy = ', airefi
	240	write(lunout,*)' nat_surf = ', nat_surf
	241	write(lunout,*)' tsurf = ', tsurf
	242	write(lunout,*)' psurf = ', psurf
	243	write(lunout,*)' zsurf = ', zsurf
	244	write(lunout,*)' rugos = ', rugos
	245	write(lunout,*)' wtsurf = ', wtsurf
	246	write(lunout,*)' wqsurf = ', wqsurf
	247	write(lunout,*)' albedo = ', albedo
	248	write(lunout,*)' xagesno = ', xagesno
	249	write(lunout,*)' restart_runoff = ', restart_runoff
	250	write(lunout,*)' qsolinp = ', qsolinp
	251	write(lunout,*)' zpicinp = ', zpicinp
	252	IF (forcing_type .eq.40) THEN
	253	write(lunout,*) '--- Forcing type GCSS Old --- with:'
	254	write(lunout,*)'imp_fcg',imp_fcg_gcssold
	255	write(lunout,*)'ts_fcg',ts_fcg_gcssold
	256	write(lunout,*)'tp_fcg',Tp_fcg_gcssold
	257	write(lunout,*)'tp_ini',Tp_ini_gcssold
	258	write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
	259	ENDIF
	260
	261	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
	262	write(lunout,*)
	263	c
	264	RETURN
	265	END
	266	!
	267	! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
	268	!
	269	c
	270	SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,
	271	& ucov,vcov,temp,q,omega2)
	272	USE dimphy
	273	USE mod_grid_phy_lmdz
	274	USE mod_phys_lmdz_para
	275	USE iophy
	276	USE phys_state_var_mod
	277	USE iostart
	278	USE write_field_phy
	279	USE infotrac
	280	use control_mod
	281
	282	IMPLICIT NONE
	283	c=======================================================
	284	c Ecriture du fichier de redemarrage sous format NetCDF
	285	c=======================================================
	286	c Declarations:
	287	c -------------
	288	#include "dimensions.h"
	289	#include "comconst.h"
	290	#include "temps.h"
	291	!!#include "control.h"
	292	#include "logic.h"
	293	#include "netcdf.inc"
	294
	295	c Arguments:
	296	c ----------
	297	CHARACTER() fichnom
	298	cAl1 plev tronque pour .nc mais plev(klev+1):=0
	299	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
	300	real :: presnivs(klon,klev)
	301	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
	302	real :: q(klon,klev,nqtot),omega2(klon,klev)
	303	c real :: ug(klev),vg(klev),fcoriolis
	304	real :: phis(klon)
	305
	306	c Variables locales pour NetCDF:
	307	c ------------------------------
	308	INTEGER nid, nvarid
	309	INTEGER idim_s
	310	INTEGER ierr, ierr_file
	311	INTEGER iq
	312	INTEGER length
	313	PARAMETER (length = 100)
	314	REAL tab_cntrl(length) ! tableau des parametres du run
	315	character*4 nmq(nqtot)
	316	character*12 modname
	317	character*80 abort_message
	318	LOGICAL found
	319	c
	320	INTEGER nb
	321
	322	modname = 'dyn1deta0 : '
	323	nmq(1)="vap"
	324	nmq(2)="cond"
	325	do iq=3,nqtot
	326	write(nmq(iq),'("tra",i1)') iq-2
	327	enddo
	328	print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
	329	CALL open_startphy(fichnom)
	330	print*,'after open startphy ',fichnom,nmq
	331
	332	c
	333	c Lecture des parametres de controle:
	334	c
	335	CALL get_var("controle",tab_cntrl)
	336
	337
	338	im = tab_cntrl(1)
	339	jm = tab_cntrl(2)
	340	lllm = tab_cntrl(3)
	341	day_ref = tab_cntrl(4)
	342	annee_ref = tab_cntrl(5)
	343	c rad = tab_cntrl(6)
	344	c omeg = tab_cntrl(7)
	345	c g = tab_cntrl(8)
	346	c cpp = tab_cntrl(9)
	347	c kappa = tab_cntrl(10)
	348	c daysec = tab_cntrl(11)
	349	c dtvr = tab_cntrl(12)
	350	c etot0 = tab_cntrl(13)
	351	c ptot0 = tab_cntrl(14)
	352	c ztot0 = tab_cntrl(15)
	353	c stot0 = tab_cntrl(16)
	354	c ang0 = tab_cntrl(17)
	355	c pa = tab_cntrl(18)
	356	c preff = tab_cntrl(19)
	357	c
	358	c clon = tab_cntrl(20)
	359	c clat = tab_cntrl(21)
	360	c grossismx = tab_cntrl(22)
	361	c grossismy = tab_cntrl(23)
	362	c
	363	IF ( tab_cntrl(24).EQ.1. ) THEN
	364	fxyhypb = . TRUE .
	365	c dzoomx = tab_cntrl(25)
	366	c dzoomy = tab_cntrl(26)
	367	c taux = tab_cntrl(28)
	368	c tauy = tab_cntrl(29)
	369	ELSE
	370	fxyhypb = . FALSE .
	371	ysinus = . FALSE .
	372	IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE.
	373	ENDIF
	374
	375	day_ini = tab_cntrl(30)
	376	itau_dyn = tab_cntrl(31)
	377	c .................................................................
	378	c
	379	c
	380	c PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
	381	cAl1
	382	Print*,'day_ref,annee_ref,day_ini,itau_dyn',
	383	& day_ref,annee_ref,day_ini,itau_dyn
	384
	385	c Lecture des champs
	386	c
	387	plev(1,klev+1)=0.
	388	CALL get_field("plev",plev,found)
	389	IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
	390	CALL get_field("play",play,found)
	391	IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
	392	CALL get_field("phi",phi,found)
	393	IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
	394	CALL get_field("phis",phis,found)
	395	IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
	396	CALL get_field("presnivs",presnivs,found)
	397	IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
	398	CALL get_field("ucov",ucov,found)
	399	IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
	400	CALL get_field("vcov",vcov,found)
	401	IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
	402	CALL get_field("temp",temp,found)
	403	IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
	404	CALL get_field("omega2",omega2,found)
	405	IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
	406
	407	Do iq=1,nqtot
	408	CALL get_field("q"//nmq(iq),q(:,:,iq),found)
	409	IF (.NOT. found)
	410	& PRINT*, modname//'Le champ <q'//nmq//'> est absent'
	411	EndDo
	412
	413	CALL close_startphy
	414	print*,' close startphy'
	415	. ,fichnom,play(1,1),play(1,klev),temp(1,klev)
	416	c
	417	RETURN
	418	END
	419	!
	420	! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
	421	!
	422	c
	423	SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,
	424	& ucov,vcov,temp,q,omega2)
	425	USE dimphy
	426	USE mod_grid_phy_lmdz
	427	USE mod_phys_lmdz_para
	428	USE phys_state_var_mod
	429	USE iostart
	430	USE infotrac
	431	use control_mod
	432
	433	IMPLICIT NONE
	434	c=======================================================
	435	c Ecriture du fichier de redemarrage sous format NetCDF
	436	c=======================================================
	437	c Declarations:
	438	c -------------
	439	#include "dimensions.h"
	440	#include "comconst.h"
	441	#include "temps.h"
	442	!!#include "control.h"
	443	#include "logic.h"
	444	#include "netcdf.inc"
	445
	446	c Arguments:
	447	c ----------
	448	CHARACTER() fichnom
	449	REAL time
	450	cAl1 plev tronque pour .nc mais plev(klev+1):=0
	451	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
	452	real :: presnivs(klon,klev)
	453	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
	454	real :: q(klon,klev,nqtot)
	455	real :: omega2(klon,klev),rho(klon,klev+1)
	456	c real :: ug(klev),vg(klev),fcoriolis
	457	real :: phis(klon)
	458
	459	c Variables locales pour NetCDF:
	460	c ------------------------------
	461	INTEGER nid, nvarid
	462	INTEGER idim_s
	463	INTEGER ierr, ierr_file
	464	INTEGER iq,l
	465	INTEGER length
	466	PARAMETER (length = 100)
	467	REAL tab_cntrl(length) ! tableau des parametres du run
	468	character*4 nmq(nqtot)
	469	character*20 modname
	470	character*80 abort_message
	471	c
	472	INTEGER nb
	473	SAVE nb
	474	DATA nb / 0 /
	475
	476	REAL zan0,zjulian,hours
	477	INTEGER yyears0,jjour0, mmois0
	478	character*30 unites
	479
	480	cDbg
	481	CALL open_restartphy(fichnom)
	482	print*,'redm1 ',fichnom,klon,klev,nqtot
	483	nmq(1)="vap"
	484	nmq(2)="cond"
	485	nmq(3)="tra1"
	486	nmq(4)="tra2"
	487
	488	modname = 'dyn1dredem'
	489	ierr = NF_OPEN(fichnom, NF_WRITE, nid)
	490	IF (ierr .NE. NF_NOERR) THEN
	491	PRINT*, "Pb. d ouverture "//fichnom
	492	CALL abort
	493	ENDIF
	494
	495	DO l=1,length
	496	tab_cntrl(l) = 0.
	497	ENDDO
	498	tab_cntrl(1) = FLOAT(iim)
	499	tab_cntrl(2) = FLOAT(jjm)
	500	tab_cntrl(3) = FLOAT(llm)
	501	tab_cntrl(4) = FLOAT(day_ref)
	502	tab_cntrl(5) = FLOAT(annee_ref)
	503	tab_cntrl(6) = rad
	504	tab_cntrl(7) = omeg
	505	tab_cntrl(8) = g
	506	tab_cntrl(9) = cpp
	507	tab_cntrl(10) = kappa
	508	tab_cntrl(11) = daysec
	509	tab_cntrl(12) = dtvr
	510	c tab_cntrl(13) = etot0
	511	c tab_cntrl(14) = ptot0
	512	c tab_cntrl(15) = ztot0
	513	c tab_cntrl(16) = stot0
	514	c tab_cntrl(17) = ang0
	515	c tab_cntrl(18) = pa
	516	c tab_cntrl(19) = preff
	517	c
	518	c ..... parametres pour le zoom ......
	519
	520	c tab_cntrl(20) = clon
	521	c tab_cntrl(21) = clat
	522	c tab_cntrl(22) = grossismx
	523	c tab_cntrl(23) = grossismy
	524	c
	525	IF ( fxyhypb ) THEN
	526	tab_cntrl(24) = 1.
	527	c tab_cntrl(25) = dzoomx
	528	c tab_cntrl(26) = dzoomy
	529	tab_cntrl(27) = 0.
	530	c tab_cntrl(28) = taux
	531	c tab_cntrl(29) = tauy
	532	ELSE
	533	tab_cntrl(24) = 0.
	534	c tab_cntrl(25) = dzoomx
	535	c tab_cntrl(26) = dzoomy
	536	tab_cntrl(27) = 0.
	537	tab_cntrl(28) = 0.
	538	tab_cntrl(29) = 0.
	539	IF( ysinus ) tab_cntrl(27) = 1.
	540	ENDIF
	541	CAl1 iday_end -> day_end
	542	tab_cntrl(30) = FLOAT(day_end)
	543	tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
	544	c
	545	CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
	546	c
	547
	548	c Ecriture/extension de la coordonnee temps
	549
	550	nb = nb + 1
	551
	552	c Ecriture des champs
	553	c
	554	CALL put_field("plev","p interfaces sauf la nulle",plev)
	555	CALL put_field("play","",play)
	556	CALL put_field("phi","geopotentielle",phi)
	557	CALL put_field("phis","geopotentiell de surface",phis)
	558	CALL put_field("presnivs","",presnivs)
	559	CALL put_field("ucov","",ucov)
	560	CALL put_field("vcov","",vcov)
	561	CALL put_field("temp","",temp)
	562	CALL put_field("omega2","",omega2)
	563
	564	Do iq=1,nqtot
	565	CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",
	566	. q(:,:,iq))
	567	EndDo
	568	CALL close_restartphy
	569
	570	c
	571	RETURN
	572	END
	573	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
	574	IMPLICIT NONE
	575	!=======================================================================
	576	! passage d'un champ de la grille scalaire a la grille physique
	577	!=======================================================================
	578
	579	!-----------------------------------------------------------------------
	580	! declarations:
	581	! -------------
	582
	583	INTEGER im,jm,ngrid,nfield
	584	REAL pdyn(im,jm,nfield)
	585	REAL pfi(ngrid,nfield)
	586
	587	INTEGER i,j,ifield,ig
	588
	589	!-----------------------------------------------------------------------
	590	! calcul:
	591	! -------
	592
	593	DO ifield=1,nfield
	594	! traitement des poles
	595	DO i=1,im
	596	pdyn(i,1,ifield)=pfi(1,ifield)
	597	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
	598	ENDDO
	599
	600	! traitement des point normaux
	601	DO j=2,jm-1
	602	ig=2+(j-2)*(im-1)
	603	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
	604	pdyn(im,j,ifield)=pdyn(1,j,ifield)
	605	ENDDO
	606	ENDDO
	607
	608	RETURN
	609	END
	610
	611
	612
	613	SUBROUTINE abort_gcm(modname, message, ierr)
	614
	615	USE IOIPSL
	616	!
	617	! Stops the simulation cleanly, closing files and printing various
	618	! comments
	619	!
	620	! Input: modname = name of calling program
	621	! message = stuff to print
	622	! ierr = severity of situation ( = 0 normal )
	623
	624	character*20 modname
	625	integer ierr
	626	character*80 message
	627
	628	write(,) 'in abort_gcm'
	629	call histclo
	630	! call histclo(2)
	631	! call histclo(3)
	632	! call histclo(4)
	633	! call histclo(5)
	634	write(,) 'out of histclo'
	635	write(,) 'Stopping in ', modname
	636	write(,) 'Reason = ',message
	637	call getin_dump
	638	!
	639	if (ierr .eq. 0) then
	640	write(,) 'Everything is cool'
	641	else
	642	write(,) 'Houston, we have a problem ', ierr
	643	endif
	644	STOP
	645	END
	646	REAL FUNCTION fq_sat(kelvin, millibar)
	647	!
	648	IMPLICIT none
	649	!======================================================================
	650	! Autheur(s): Z.X. Li (LMD/CNRS)
	651	! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
	652	!======================================================================
	653	! Arguments:
	654	! kelvin---input-R: temperature en Kelvin
	655	! millibar--input-R: pression en mb
	656	!
	657	! fq_sat----output-R: vapeur d'eau saturante en kg/kg
	658	!======================================================================
	659	!
	660	REAL kelvin, millibar
	661	!
	662	REAL r2es
	663	PARAMETER (r2es=611.14 *18.0153/28.9644)
	664	!
	665	REAL r3les, r3ies, r3es
	666	PARAMETER (R3LES=17.269)
	667	PARAMETER (R3IES=21.875)
	668	!
	669	REAL r4les, r4ies, r4es
	670	PARAMETER (R4LES=35.86)
	671	PARAMETER (R4IES=7.66)
	672	!
	673	REAL rtt
	674	PARAMETER (rtt=273.16)
	675	!
	676	REAL retv
	677	PARAMETER (retv=28.9644/18.0153 - 1.0)
	678	!
	679	REAL zqsat
	680	REAL temp, pres
	681	! ------------------------------------------------------------------
	682	!
	683	!
	684	temp = kelvin
	685	pres = millibar * 100.0
	686	! write(,)'kelvin,millibar=',kelvin,millibar
	687	! write(,)'temp,pres=',temp,pres
	688	!
	689	IF (temp .LE. rtt) THEN
	690	r3es = r3ies
	691	r4es = r4ies
	692	ELSE
	693	r3es = r3les
	694	r4es = r4les
	695	ENDIF
	696	!
	697	zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
	698	zqsat=MIN(0.5,ZQSAT)
	699	zqsat=zqsat/(1.-retv *zqsat)
	700	!
	701	fq_sat = zqsat
	702	!
	703	RETURN
	704	END
	705	subroutine scopy(n,sx,incx,sy,incy)
	706	!
	707	IMPLICIT NONE
	708	!
	709	integer n,incx,incy,ix,iy,i
	710	real sx((n-1)incx+1),sy((n-1)incy+1)
	711	!
	712	iy=1
	713	ix=1
	714	do 10 i=1,n
	715	sy(iy)=sx(ix)
	716	ix=ix+incx
	717	iy=iy+incy
	718	10 continue
	719	!
	720	return
	721	end
	722	subroutine wrgradsfi(if,nl,field,name,titlevar)
	723	implicit none
	724
	725	! Declarations
	726
	727	#include "gradsdef.h"
	728
	729	! arguments
	730	integer if,nl
	731	real field(imxjmxlmx)
	732	character*10 name,file
	733	character*10 titlevar
	734
	735	! local
	736
[1608]	737	integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
[1607]	738
	739	logical writectl
	740
	741
	742	writectl=.false.
	743
	744	! print*,if,iid(if),jid(if),ifd(if),jfd(if)
	745	iii=iid(if)
	746	iji=jid(if)
	747	iif=ifd(if)
	748	ijf=jfd(if)
	749	im=iif-iii+1
	750	jm=ijf-iji+1
	751	lm=lmd(if)
	752
	753
	754	! print*,'im,jm,lm,name,firsttime(if)'
	755	! print*,im,jm,lm,name,firsttime(if)
	756
	757	if(firsttime(if)) then
	758	if(name.eq.var(1,if)) then
	759	firsttime(if)=.false.
	760	ivar(if)=1
	761	print*,'fin de l initialiation de l ecriture du fichier'
	762	print*,file
	763	print*,'fichier no: ',if
	764	print*,'unit ',unit(if)
	765	print*,'nvar ',nvar(if)
	766	print*,'vars ',(var(iv,if),iv=1,nvar(if))
	767	else
	768	ivar(if)=ivar(if)+1
	769	nvar(if)=ivar(if)
	770	var(ivar(if),if)=name
[1608]	771	tvar(ivar(if),if)=trim(titlevar)
[1607]	772	nld(ivar(if),if)=nl
	773	print*,'initialisation ecriture de ',var(ivar(if),if)
	774	print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
	775	endif
	776	writectl=.true.
	777	itime(if)=1
	778	else
	779	ivar(if)=mod(ivar(if),nvar(if))+1
	780	if (ivar(if).eq.nvar(if)) then
	781	writectl=.true.
	782	itime(if)=itime(if)+1
	783	endif
	784
	785	if(var(ivar(if),if).ne.name) then
	786	print*,'Il faut stoker la meme succession de champs a chaque'
	787	print*,'pas de temps'
	788	print*,'fichier no: ',if
	789	print*,'unit ',unit(if)
	790	print*,'nvar ',nvar(if)
	791	print*,'vars ',(var(iv,if),iv=1,nvar(if))
	792
	793	stop
	794	endif
	795	endif
	796
	797	! print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
	798	! print*,ivar(if),nvar(if),var(ivar(if),if),writectl
	799	do l=1,nl
	800	irec(if)=irec(if)+1
	801	! print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf,
	802	! s (l-1)imd(if)jmd(if)+(iji-1)*imd(if)+iii
	803	! s ,(l-1)imd(if)jmd(if)+(ijf-1)*imd(if)+iif
	804	write(unit(if)+1,rec=irec(if))
	805	s ((field((l-1)imd(if)jmd(if)+(j-1)*imd(if)+i)
	806	s ,i=iii,iif),j=iji,ijf)
	807	enddo
	808	if (writectl) then
	809
	810	file=fichier(if)
	811	! WARNING! on reecrase le fichier .ctl a chaque ecriture
[1608]	812	open(unit(if),file=trim(file)//'.ctl',
[1607]	813	& form='formatted',status='unknown')
	814	write(unit(if),'(a5,1x,a40)')
[1608]	815	& 'DSET ','^'//trim(file)//'.dat'
[1607]	816
	817	write(unit(if),'(a12)') 'UNDEF 1.0E30'
	818	write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if)
	819	call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF')
	820	call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF')
	821	call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF')
	822	write(unit(if),'(a4,i10,a30)')
	823	& 'TDEF ',itime(if),' LINEAR 07AUG1998 30MN '
	824	write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if)
	825	do iv=1,nvar(if)
	826	! print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)'
	827	! print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)
	828	write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if)
	829	& ,99,tvar(iv,if)
	830	enddo
	831	write(unit(if),'(a7)') 'ENDVARS'
	832	!
	833	1000 format(a5,3x,i4,i3,1x,a39)
	834
	835	close(unit(if))
	836
	837	endif ! writectl
	838
	839	return
	840
	841	END
	842
	843	subroutine inigrads(if,im
	844	s ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz
	845	s ,dt,file,titlel)
	846
	847
	848	implicit none
	849
[1608]	850	integer if,im,jm,lm,i,j,l
[1607]	851	real x(im),y(jm),z(lm),fx,fy,fz,dt
	852	real xmin,xmax,ymin,ymax
	853	integer nf
	854
	855	character file10,titlel40
	856
	857	#include "gradsdef.h"
	858
	859	data unit/24,32,34,36,38,40,42,44,46,48/
	860	data nf/0/
	861
	862	if (if.le.nf) stop'verifier les appels a inigrads'
	863
	864	print*,'Entree dans inigrads'
	865
	866	nf=if
	867	title(if)=titlel
	868	ivar(if)=0
	869
[1608]	870	fichier(if)=trim(file)
[1607]	871
	872	firsttime(if)=.true.
	873	dtime(if)=dt
	874
	875	iid(if)=1
	876	ifd(if)=im
	877	imd(if)=im
	878	do i=1,im
	879	xd(i,if)=x(i)*fx
	880	if(xd(i,if).lt.xmin) iid(if)=i+1
	881	if(xd(i,if).le.xmax) ifd(if)=i
	882	enddo
	883	print*,'On stoke du point ',iid(if),' a ',ifd(if),' en x'
	884
	885	jid(if)=1
	886	jfd(if)=jm
	887	jmd(if)=jm
	888	do j=1,jm
	889	yd(j,if)=y(j)*fy
	890	if(yd(j,if).gt.ymax) jid(if)=j+1
	891	if(yd(j,if).ge.ymin) jfd(if)=j
	892	enddo
	893	print*,'On stoke du point ',jid(if),' a ',jfd(if),' en y'
	894
	895	print*,'Open de dat'
	896	print*,'file=',file
	897	print*,'fichier(if)=',fichier(if)
	898
	899	print,4(ifd(if)-iid(if))*(jfd(if)-jid(if))
[1608]	900	print*,trim(file)//'.dat'
[1607]	901
[1608]	902	OPEN (unit(if)+1,FILE=trim(file)//'.dat',
[1607]	903	s FORM='UNFORMATTED',
	904	s ACCESS='DIRECT'
	905	s ,RECL=4(ifd(if)-iid(if)+1)(jfd(if)-jid(if)+1))
	906
	907	print*,'Open de dat ok'
	908
	909	lmd(if)=lm
	910	do l=1,lm
	911	zd(l,if)=z(l)*fz
	912	enddo
	913
	914	irec(if)=0
	915	!CR
	916	! print*,if,imd(if),jmd(if),lmd(if)
	917	! print*,'if,imd(if),jmd(if),lmd(if)'
	918
	919	return
	920	end
	921	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
	922	IMPLICIT NONE
	923	!=======================================================================
	924	! passage d'un champ de la grille scalaire a la grille physique
	925	!=======================================================================
	926
	927	!-----------------------------------------------------------------------
	928	! declarations:
	929	! -------------
	930
	931	INTEGER im,jm,ngrid,nfield
	932	REAL pdyn(im,jm,nfield)
	933	REAL pfi(ngrid,nfield)
	934
	935	INTEGER j,ifield,ig
	936
	937	!-----------------------------------------------------------------------
	938	! calcul:
	939	! -------
	940
	941	IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
	942	s STOP 'probleme de dim'
	943	! traitement des poles
	944	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
	945	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
	946
	947	! traitement des point normaux
	948	DO ifield=1,nfield
	949	DO j=2,jm-1
	950	ig=2+(j-2)*(im-1)
	951	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
	952	ENDDO
	953	ENDDO
	954
	955	RETURN
	956	END
	957
	958	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	959	subroutine writelim
	960	s (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
	961	s phy_fter,phy_foce,phy_flic,phy_fsic)
	962
	963	use dimphy
	964	implicit none
	965	!
	966	#include "dimensions.h"
	967	#include "netcdf.inc"
	968
	969	integer klon1d
	970	integer k
	971	parameter (klon1d=1)
	972	REAL phy_nat(klon1d,360)
	973	REAL phy_alb(klon1d,360)
	974	REAL phy_sst(klon1d,360)
	975	REAL phy_bil(klon1d,360)
	976	REAL phy_rug(klon1d,360)
	977	REAL phy_ice(klon1d,360)
	978	REAL phy_fter(klon1d,360)
	979	REAL phy_foce(klon1d,360)
	980	REAL phy_flic(klon1d,360)
	981	REAL phy_fsic(klon1d,360)
	982
	983	INTEGER ierr
	984	INTEGER dimfirst(3)
	985	INTEGER dimlast(3)
	986	!
	987	INTEGER nid, ndim, ntim
	988	INTEGER dims(2), debut(2), epais(2)
	989	INTEGER id_tim
	990	INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
	991	INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
	992
	993	PRINT*, 'Ecriture du fichier limit'
	994	!
	995	ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
	996	!
	997	ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
	998	. "Fichier conditions aux limites")
	999	ierr = NF_DEF_DIM (nid, "points_physiques", klon1d, ndim)
	1000	ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
	1001	!
	1002	dims(1) = ndim
	1003	dims(2) = ntim
	1004	!
	1005	!cc ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
	1006	ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
	1007	ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
	1008	. "Jour dans l annee")
	1009	!cc ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
	1010	ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
	1011	ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
	1012	. "Nature du sol (0,1,2,3)")
	1013	!cc ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
	1014	ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
	1015	ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
	1016	. "Temperature superficielle de la mer")
	1017	!cc ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
	1018	ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
	1019	ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
	1020	. "Reference flux de chaleur au sol")
	1021	!cc ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
	1022	ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
	1023	ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
	1024	. "Albedo a la surface")
	1025	!cc ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
	1026	ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
	1027	ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
	1028	. "Rugosite")
	1029
	1030	ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
	1031	ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
	1032	ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
	1033	ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
	1034	ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
	1035	ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
	1036	ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
	1037	ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
	1038	!
	1039	ierr = NF_ENDDEF(nid)
	1040	!
	1041	DO k = 1, 360
	1042	!
	1043	debut(1) = 1
	1044	debut(2) = k
	1045	epais(1) = klon1d
	1046	epais(2) = 1
	1047	!
	1048	#ifdef NC_DOUBLE
	1049	ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
	1050	ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
	1051	ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
	1052	ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
	1053	ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
	1054	ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
	1055	ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k))
	1056	ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k))
	1057	ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k))
	1058	ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k))
	1059	#else
	1060	ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
	1061	ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
	1062	ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
	1063	ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
	1064	ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
	1065	ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
	1066	ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k))
	1067	ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k))
	1068	ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k))
	1069	ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k))
	1070
	1071	#endif
	1072	!
	1073	ENDDO
	1074	!
	1075	ierr = NF_CLOSE(nid)
	1076	!
	1077	return
	1078	end
	1079
	1080	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1081	SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
	1082
	1083	! Auteur : P. Le Van .
	1084	!
	1085	IMPLICIT NONE
	1086
	1087	#include "dimensions.h"
	1088	#include "paramet.h"
	1089	!
	1090	!=======================================================================
	1091	!
	1092	!
	1093	! s = sigma ** kappa : coordonnee verticale
	1094	! dsig(l) : epaisseur de la couche l ds la coord. s
	1095	! sig(l) : sigma a l'interface des couches l et l-1
	1096	! ds(l) : distance entre les couches l et l-1 en coord.s
	1097	!
	1098	!=======================================================================
	1099	!
	1100	REAL pa,preff
	1101	REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
	1102	REAL presnivs(llm)
	1103	!
	1104	! declarations:
	1105	! -------------
	1106	!
	1107	REAL sig(llm+1),dsig(llm)
	1108	!
	1109	INTEGER l
	1110	REAL snorm
	1111	REAL alpha,beta,gama,delta,deltaz,h
	1112	INTEGER np,ierr
	1113	REAL pi,x
	1114
	1115	!-----------------------------------------------------------------------
	1116	!
	1117	pi=2.*ASIN(1.)
	1118
	1119	OPEN(99,file='sigma.def',status='old',form='formatted',
	1120	s iostat=ierr)
	1121
	1122	!-----------------------------------------------------------------------
	1123	! cas 1 on lit les options dans sigma.def:
	1124	! ----------------------------------------
	1125
	1126	IF (ierr.eq.0) THEN
	1127
	1128	print*,'WARNING!!! on lit les options dans sigma.def'
	1129	READ(99,*) deltaz
	1130	READ(99,*) h
	1131	READ(99,*) beta
	1132	READ(99,*) gama
	1133	READ(99,*) delta
	1134	READ(99,*) np
	1135	CLOSE(99)
	1136	alpha=deltaz/(llm*h)
	1137	!
	1138
	1139	DO 1 l = 1, llm
	1140	dsig(l) = (alpha+(1.-alpha)exp(-beta(llm-l)))*
	1141	$ ( (tanh(gamal)/tanh(gamallm))**np +
	1142	$ (1.-l/FLOAT(llm))*delta )
	1143	1 CONTINUE
	1144
	1145	sig(1)=1.
	1146	DO 101 l=1,llm-1
	1147	sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
	1148	101 CONTINUE
	1149	sig(llm+1)=0.
	1150
	1151	DO 2 l = 1, llm
	1152	dsig(l) = sig(l)-sig(l+1)
	1153	2 CONTINUE
	1154	!
	1155
	1156	ELSE
	1157	!-----------------------------------------------------------------------
	1158	! cas 2 ancienne discretisation (LMD5...):
	1159	! ----------------------------------------
	1160
	1161	PRINT*,'WARNING!!! Ancienne discretisation verticale'
	1162
	1163	h=7.
	1164	snorm = 0.
	1165	DO l = 1, llm
	1166	x = 2.asin(1.) (FLOAT(l)-0.5) / float(llm+1)
	1167	dsig(l) = 1.0 + 7.0 * SIN(x)**2
	1168	snorm = snorm + dsig(l)
	1169	ENDDO
	1170	snorm = 1./snorm
	1171	DO l = 1, llm
	1172	dsig(l) = dsig(l)*snorm
	1173	ENDDO
	1174	sig(llm+1) = 0.
	1175	DO l = llm, 1, -1
	1176	sig(l) = sig(l+1) + dsig(l)
	1177	ENDDO
	1178
	1179	ENDIF
	1180
	1181
	1182	DO l=1,llm
	1183	nivsigs(l) = FLOAT(l)
	1184	ENDDO
	1185
	1186	DO l=1,llmp1
	1187	nivsig(l)= FLOAT(l)
	1188	ENDDO
	1189
	1190	!
	1191	! .... Calculs de ap(l) et de bp(l) ....
	1192	! .........................................
	1193	!
	1194	!
	1195	! ..... pa et preff sont lus sur les fichiers start par lectba .....
	1196	!
	1197
	1198	bp(llmp1) = 0.
	1199
	1200	DO l = 1, llm
	1201	!c
	1202	!cc ap(l) = 0.
	1203	!cc bp(l) = sig(l)
	1204
	1205	bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
	1206	ap(l) = pa * ( sig(l) - bp(l) )
	1207	!
	1208	ENDDO
	1209	ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
	1210
	1211	PRINT *,' BP '
	1212	PRINT *, bp
	1213	PRINT *,' AP '
	1214	PRINT *, ap
	1215
	1216	DO l = 1, llm
	1217	dpres(l) = bp(l) - bp(l+1)
	1218	presnivs(l) = 0.5 ( ap(l)+bp(l)preff + ap(l+1)+bp(l+1)*preff )
	1219	ENDDO
	1220
	1221	PRINT *,' PRESNIVS '
	1222	PRINT *,presnivs
	1223
	1224	RETURN
	1225	END
	1226
	1227	!======================================================================
	1228	SUBROUTINE read_tsurf1d(knon,knindex,sst_out)
	1229
	1230	! This subroutine specifies the surface temperature to be used in 1D simulations
	1231
	1232	USE dimphy, ONLY : klon
	1233
	1234	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
	1235	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
	1236	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
	1237
	1238	INTEGER :: i
	1239	! COMMON defined in lmdz1d.F:
	1240	real ts_cur
	1241	common /sst_forcing/ts_cur
	1242
	1243	DO i = 1, knon
	1244	sst_out(i) = ts_cur
	1245	ENDDO
	1246
	1247	END SUBROUTINE read_tsurf1d
	1248
	1249	!===============================================================
	1250	subroutine advect_vert(llm,w,dt,q,plev)
	1251	!===============================================================
	1252	! Schema amont pour l'advection verticale en 1D
	1253	! w est la vitesse verticale dp/dt en Pa/s
	1254	! Traitement en volumes finis
	1255	! d / dt ( zm q ) = delta_z ( omega q )
	1256	! d / dt ( zm ) = delta_z ( omega )
	1257	! avec zm = delta_z ( p )
	1258	! si * designe la valeur au pas de temps t+dt
	1259	! zm(l) q(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
	1260	! zm*(l) -zm(l) = w(l+1) - w(l)
	1261	! avec w=omega * dt
	1262	!---------------------------------------------------------------
	1263	implicit none
	1264	! arguments
	1265	integer llm
	1266	real w(llm+1),q(llm),plev(llm+1),dt
	1267
	1268	! local
	1269	integer l
	1270	real zwq(llm+1),zm(llm+1),zw(llm+1)
	1271	real qold
	1272
	1273	!---------------------------------------------------------------
	1274
	1275	do l=1,llm
	1276	zw(l)=dt*w(l)
	1277	zm(l)=plev(l)-plev(l+1)
	1278	zwq(l)=q(l)*zw(l)
	1279	enddo
	1280	zwq(llm+1)=0.
	1281	zw(llm+1)=0.
	1282
	1283	do l=1,llm
	1284	qold=q(l)
	1285	q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
	1286	print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
	1287	enddo
	1288
	1289
	1290	return
	1291	end
	1292
	1293	!===============================================================
	1294
	1295
	1296	SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
	1297	! q,temp,u,v,
	1298	! play,plev)
	1299	!itlmd
	1300	!----------------------------------------------------------------------
	1301	! Calcul de l'advection verticale (ascendance et subsidence) de
	1302	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
	1303	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
	1304	! sans WTG rajouter une advection horizontale
	1305	!----------------------------------------------------------------------
	1306	implicit none
	1307	#include "YOMCST.h"
	1308	! argument
	1309	integer llm
	1310	real omega(llm+1),d_t_va(llm), d_q_va(llm,3)
	1311	real d_u_va(llm), d_v_va(llm)
	1312	real q(llm,3),temp(llm)
	1313	real u(llm),v(llm)
	1314	real play(llm),plev(llm+1)
	1315	! interne
	1316	integer l
	1317	real alpha,omgdown,omgup
	1318
	1319	do l= 1,llm
	1320	if(l.eq.1) then
	1321	!si omgup pour la couche 1, alors tendance nulle
	1322	omgdown=max(omega(2),0.0)
	1323	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
	1324	& (1.+q(l,1)))
	1325	d_t_va(l)= alpha*(omgdown)-
	1326	& omgdown*(temp(l)-temp(l+1))
	1327	& /(play(l)-play(l+1))
	1328
	1329	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
	1330	& /(play(l)-play(l+1))
	1331
	1332	d_u_va(l)= -omgdown*(u(l)-u(l+1))
	1333	& /(play(l)-play(l+1))
	1334	d_v_va(l)= -omgdown*(v(l)-v(l+1))
	1335	& /(play(l)-play(l+1))
	1336
	1337
	1338	elseif(l.eq.llm) then
	1339	omgup=min(omega(l),0.0)
	1340	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
	1341	& (1.+q(l,1)))
	1342	d_t_va(l)= alpha*(omgup)-
	1343	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
	1344	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
	1345	d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
	1346	d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
	1347	d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
	1348
	1349	else
	1350	omgup=min(omega(l),0.0)
	1351	omgdown=max(omega(l+1),0.0)
	1352	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
	1353	& (1.+q(l,1)))
	1354	d_t_va(l)= alpha*(omgup+omgdown)-
	1355	& omgdown*(temp(l)-temp(l+1))
	1356	& /(play(l)-play(l+1))-
	1357	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
	1358	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
	1359	! print*, ' ??? '
	1360
	1361	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
	1362	& /(play(l)-play(l+1))-
	1363	& omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
	1364	d_u_va(l)= -omgdown*(u(l)-u(l+1))
	1365	& /(play(l)-play(l+1))-
	1366	& omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
	1367	d_v_va(l)= -omgdown*(v(l)-v(l+1))
	1368	& /(play(l)-play(l+1))-
	1369	& omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
	1370
	1371	endif
	1372
	1373	enddo
	1374	!fin itlmd
	1375	return
	1376	end
	1377	! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
	1378	SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,
	1379	! q,temp,u,v,play)
	1380	!itlmd
	1381	!----------------------------------------------------------------------
	1382	! Calcul de l'advection verticale (ascendance et subsidence) de
	1383	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
	1384	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
	1385	! sans WTG rajouter une advection horizontale
	1386	!----------------------------------------------------------------------
	1387	implicit none
	1388	#include "YOMCST.h"
	1389	! argument
	1390	integer llm,nqtot
	1391	real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
	1392	! real d_u_va(llm), d_v_va(llm)
	1393	real q(llm,nqtot),temp(llm)
	1394	real u(llm),v(llm)
	1395	real play(llm)
	1396	real cor(llm)
	1397	! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
	1398	real dph(llm),dqdp(llm),dtdp(llm)
	1399	! interne
	1400	integer l,k
	1401	real alpha,omdn,omup
	1402
	1403	! dudp=0.
	1404	! dvdp=0.
	1405	dqdp=0.
	1406	dtdp=0.
	1407	! d_u_va=0.
	1408	! d_v_va=0.
	1409
	1410	cor(:) = rkappatemp(1.+q(:,1)rv/rd)/(play(1.+q(:,1)))
	1411
	1412
	1413	do k=2,llm-1
	1414
	1415	dph (k-1) = (play(k )- play(k-1 ))
	1416	! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
	1417	! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1)
	1418	dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1)
	1419	dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1)
	1420
	1421	enddo
	1422
	1423	! dudp ( llm ) = dudp ( llm-1 )
	1424	! dvdp ( llm ) = dvdp ( llm-1 )
	1425	dqdp ( llm ) = dqdp ( llm-1 )
	1426	dtdp ( llm ) = dtdp ( llm-1 )
	1427
	1428	do k=2,llm-1
	1429	omdn=max(0.0,omega(k+1))
	1430	omup=min(0.0,omega( k ))
	1431
	1432	! d_u_va(k) = -omdndudp(k)-omupdudp(k-1)
	1433	! d_v_va(k) = -omdndvdp(k)-omupdvdp(k-1)
	1434	d_q_va(k,1)= -omdndqdp(k)-omupdqdp(k-1)
	1435	d_t_va(k) = -omdndtdp(k)-omupdtdp(k-1)
	1436	: +(omup+omdn)*cor(k)
	1437	enddo
	1438
	1439	omdn=max(0.0,omega( 2 ))
	1440	omup=min(0.0,omega(llm))
	1441	! d_u_va( 1 ) = -omdn*dudp( 1 )
	1442	! d_u_va(llm) = -omup*dudp(llm)
	1443	! d_v_va( 1 ) = -omdn*dvdp( 1 )
	1444	! d_v_va(llm) = -omup*dvdp(llm)
	1445	d_q_va( 1 ,1) = -omdn*dqdp( 1 )
	1446	d_q_va(llm,1) = -omup*dqdp(llm)
	1447	d_t_va( 1 ) = -omdndtdp( 1 )+omdncor( 1 )
	1448	d_t_va(llm) = -omupdtdp(llm)!+omupcor(llm)
	1449
	1450	! if(abs(rlat(1))>10.) then
	1451	! Calculate the tendency due agestrophic motions
	1452	! du_age = fcoriolis*(v-vg)
	1453	! dv_age = fcoriolis*(ug-u)
	1454	! endif
	1455
	1456	! call writefield_phy('d_t_va',d_t_va,llm)
	1457
	1458	return
	1459	end
	1460
	1461	!======================================================================
	1462	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga
	1463	: ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
	1464	: ,ht_toga,vt_toga,hq_toga,vq_toga)
	1465	implicit none
	1466
	1467	c-------------------------------------------------------------------------
	1468	c Read TOGA-COARE forcing data
	1469	c-------------------------------------------------------------------------
	1470
	1471	integer nlev_toga,nt_toga
	1472	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
	1473	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
	1474	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
	1475	real w_toga(nlev_toga,nt_toga)
	1476	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
	1477	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
	1478	character*80 fich_toga
	1479
	1480	integer no,l,k,ip
	1481	real riy,rim,rid,rih,bid
	1482
	1483	integer iy,im,id,ih
	1484
	1485	real plev_min
	1486
	1487	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
	1488
[1608]	1489	open(21,file=trim(fich_toga),form='formatted')
[1607]	1490	read(21,'(a)')
	1491	do ip = 1, nt_toga
	1492	read(21,'(a)')
	1493	read(21,'(a)')
	1494	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
	1495	read(21,'(a)')
	1496	read(21,'(a)')
	1497
	1498	do k = 1, nlev_toga
	1499	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)
	1500	: ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)
	1501	: ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
	1502
	1503	! conversion in SI units:
	1504	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
	1505	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
	1506	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
	1507	! no water vapour tendency above 55 hPa
	1508	if (plev_toga(k,ip) .lt. plev_min) then
	1509	q_toga(k,ip) = 0.
	1510	hq_toga(k,ip) = 0.
	1511	vq_toga(k,ip) =0.
	1512	endif
	1513	enddo
	1514
	1515	ts_toga(ip)=ts_toga(ip)+273.15 ! K
	1516	enddo
	1517	close(21)
	1518
	1519	223 format(4i3,6f8.2)
	1520	226 format(f7.1,1x,10f8.2)
	1521	227 format(f7.1,1x,1p,4e11.3)
	1522	230 format(6f9.3,4e11.3)
	1523
	1524	return
	1525	end
	1526
	1527	!=====================================================================
	1528	subroutine read_twpice(fich_twpice,nlevel,ntime
	1529	: ,T_srf,plev,T,q,u,v,omega
	1530	: ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
	1531
	1532	!program reading forcings of the TWP-ICE experiment
	1533
	1534	! use netcdf
	1535
	1536	implicit none
	1537
	1538	#include "netcdf.inc"
	1539
	1540	integer ntime,nlevel
	1541	integer l,k
	1542	character*80 :: fich_twpice
	1543	real*8 time(ntime)
	1544	real*8 lat, lon, alt, phis
	1545	real*8 lev(nlevel)
	1546	real*8 plev(nlevel,ntime)
	1547
	1548	real*8 T(nlevel,ntime)
	1549	real*8 q(nlevel,ntime),u(nlevel,ntime)
	1550	real*8 v(nlevel,ntime)
	1551	real*8 omega(nlevel,ntime), div(nlevel,ntime)
	1552	real*8 T_adv_h(nlevel,ntime)
	1553	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
	1554	real*8 q_adv_v(nlevel,ntime)
	1555	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
	1556	real*8 s_adv_v(nlevel,ntime)
	1557	real*8 p_srf_aver(ntime), p_srf_center(ntime)
	1558	real*8 T_srf(ntime)
	1559
	1560	integer nid, ierr
	1561	integer nbvar3d
	1562	parameter(nbvar3d=20)
	1563	integer var3didin(nbvar3d)
	1564
	1565	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
	1566	if (ierr.NE.NF_NOERR) then
	1567	write(,) 'ERROR: Pb opening forcings cdf file '
	1568	write(,) NF_STRERROR(ierr)
	1569	stop ""
	1570	endif
	1571
	1572	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
	1573	if(ierr/=NF_NOERR) then
	1574	write(,) NF_STRERROR(ierr)
	1575	stop 'lat'
	1576	endif
	1577
	1578	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
	1579	if(ierr/=NF_NOERR) then
	1580	write(,) NF_STRERROR(ierr)
	1581	stop 'lon'
	1582	endif
	1583
	1584	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
	1585	if(ierr/=NF_NOERR) then
	1586	write(,) NF_STRERROR(ierr)
	1587	stop 'alt'
	1588	endif
	1589
	1590	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
	1591	if(ierr/=NF_NOERR) then
	1592	write(,) NF_STRERROR(ierr)
	1593	stop 'phis'
	1594	endif
	1595
	1596	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
	1597	if(ierr/=NF_NOERR) then
	1598	write(,) NF_STRERROR(ierr)
	1599	stop 'T'
	1600	endif
	1601
	1602	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
	1603	if(ierr/=NF_NOERR) then
	1604	write(,) NF_STRERROR(ierr)
	1605	stop 'q'
	1606	endif
	1607
	1608	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
	1609	if(ierr/=NF_NOERR) then
	1610	write(,) NF_STRERROR(ierr)
	1611	stop 'u'
	1612	endif
	1613
	1614	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
	1615	if(ierr/=NF_NOERR) then
	1616	write(,) NF_STRERROR(ierr)
	1617	stop 'v'
	1618	endif
	1619
	1620	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
	1621	if(ierr/=NF_NOERR) then
	1622	write(,) NF_STRERROR(ierr)
	1623	stop 'omega'
	1624	endif
	1625
	1626	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
	1627	if(ierr/=NF_NOERR) then
	1628	write(,) NF_STRERROR(ierr)
	1629	stop 'div'
	1630	endif
	1631
	1632	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
	1633	if(ierr/=NF_NOERR) then
	1634	write(,) NF_STRERROR(ierr)
	1635	stop 'T_adv_h'
	1636	endif
	1637
	1638	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
	1639	if(ierr/=NF_NOERR) then
	1640	write(,) NF_STRERROR(ierr)
	1641	stop 'T_adv_v'
	1642	endif
	1643
	1644	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
	1645	if(ierr/=NF_NOERR) then
	1646	write(,) NF_STRERROR(ierr)
	1647	stop 'q_adv_h'
	1648	endif
	1649
	1650	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
	1651	if(ierr/=NF_NOERR) then
	1652	write(,) NF_STRERROR(ierr)
	1653	stop 'q_adv_v'
	1654	endif
	1655
	1656	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
	1657	if(ierr/=NF_NOERR) then
	1658	write(,) NF_STRERROR(ierr)
	1659	stop 's'
	1660	endif
	1661
	1662	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
	1663	if(ierr/=NF_NOERR) then
	1664	write(,) NF_STRERROR(ierr)
	1665	stop 's_adv_h'
	1666	endif
	1667
	1668	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
	1669	if(ierr/=NF_NOERR) then
	1670	write(,) NF_STRERROR(ierr)
	1671	stop 's_adv_v'
	1672	endif
	1673
	1674	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
	1675	if(ierr/=NF_NOERR) then
	1676	write(,) NF_STRERROR(ierr)
	1677	stop 'p_srf_aver'
	1678	endif
	1679
	1680	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
	1681	if(ierr/=NF_NOERR) then
	1682	write(,) NF_STRERROR(ierr)
	1683	stop 'p_srf_center'
	1684	endif
	1685
	1686	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
	1687	if(ierr/=NF_NOERR) then
	1688	write(,) NF_STRERROR(ierr)
	1689	stop 'T_srf'
	1690	endif
	1691
	1692	!dimensions lecture
	1693	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
	1694
	1695	!pressure
	1696	do l=1,ntime
	1697	do k=1,nlevel
	1698	plev(k,l)=lev(k)
	1699	enddo
	1700	enddo
	1701
	1702	#ifdef NC_DOUBLE
	1703	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
	1704	#else
	1705	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
	1706	#endif
	1707	if(ierr/=NF_NOERR) then
	1708	write(,) NF_STRERROR(ierr)
	1709	stop "getvarup"
	1710	endif
	1711	! write(,)'lecture lat ok',lat
	1712
	1713	#ifdef NC_DOUBLE
	1714	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
	1715	#else
	1716	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
	1717	#endif
	1718	if(ierr/=NF_NOERR) then
	1719	write(,) NF_STRERROR(ierr)
	1720	stop "getvarup"
	1721	endif
	1722	! write(,)'lecture lon ok',lon
	1723
	1724	#ifdef NC_DOUBLE
	1725	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
	1726	#else
	1727	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
	1728	#endif
	1729	if(ierr/=NF_NOERR) then
	1730	write(,) NF_STRERROR(ierr)
	1731	stop "getvarup"
	1732	endif
	1733	! write(,)'lecture alt ok',alt
	1734
	1735	#ifdef NC_DOUBLE
	1736	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
	1737	#else
	1738	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
	1739	#endif
	1740	if(ierr/=NF_NOERR) then
	1741	write(,) NF_STRERROR(ierr)
	1742	stop "getvarup"
	1743	endif
	1744	! write(,)'lecture phis ok',phis
	1745
	1746	#ifdef NC_DOUBLE
	1747	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
	1748	#else
	1749	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
	1750	#endif
	1751	if(ierr/=NF_NOERR) then
	1752	write(,) NF_STRERROR(ierr)
	1753	stop "getvarup"
	1754	endif
	1755	! write(,)'lecture T ok'
	1756
	1757	#ifdef NC_DOUBLE
	1758	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
	1759	#else
	1760	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
	1761	#endif
	1762	if(ierr/=NF_NOERR) then
	1763	write(,) NF_STRERROR(ierr)
	1764	stop "getvarup"
	1765	endif
	1766	! write(,)'lecture q ok'
	1767	!q in kg/kg
	1768	do l=1,ntime
	1769	do k=1,nlevel
	1770	q(k,l)=q(k,l)/1000.
	1771	enddo
	1772	enddo
	1773	#ifdef NC_DOUBLE
	1774	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
	1775	#else
	1776	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
	1777	#endif
	1778	if(ierr/=NF_NOERR) then
	1779	write(,) NF_STRERROR(ierr)
	1780	stop "getvarup"
	1781	endif
	1782	! write(,)'lecture u ok'
	1783
	1784	#ifdef NC_DOUBLE
	1785	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
	1786	#else
	1787	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
	1788	#endif
	1789	if(ierr/=NF_NOERR) then
	1790	write(,) NF_STRERROR(ierr)
	1791	stop "getvarup"
	1792	endif
	1793	! write(,)'lecture v ok'
	1794
	1795	#ifdef NC_DOUBLE
	1796	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
	1797	#else
	1798	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
	1799	#endif
	1800	if(ierr/=NF_NOERR) then
	1801	write(,) NF_STRERROR(ierr)
	1802	stop "getvarup"
	1803	endif
	1804	! write(,)'lecture omega ok'
	1805	!omega in mb/hour
	1806	do l=1,ntime
	1807	do k=1,nlevel
	1808	omega(k,l)=omega(k,l)*100./3600.
	1809	enddo
	1810	enddo
	1811
	1812	#ifdef NC_DOUBLE
	1813	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
	1814	#else
	1815	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
	1816	#endif
	1817	if(ierr/=NF_NOERR) then
	1818	write(,) NF_STRERROR(ierr)
	1819	stop "getvarup"
	1820	endif
	1821	! write(,)'lecture div ok'
	1822
	1823	#ifdef NC_DOUBLE
	1824	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
	1825	#else
	1826	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
	1827	#endif
	1828	if(ierr/=NF_NOERR) then
	1829	write(,) NF_STRERROR(ierr)
	1830	stop "getvarup"
	1831	endif
	1832	! write(,)'lecture T_adv_h ok'
	1833	!T adv in K/s
	1834	do l=1,ntime
	1835	do k=1,nlevel
	1836	T_adv_h(k,l)=T_adv_h(k,l)/3600.
	1837	enddo
	1838	enddo
	1839
	1840
	1841	#ifdef NC_DOUBLE
	1842	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
	1843	#else
	1844	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
	1845	#endif
	1846	if(ierr/=NF_NOERR) then
	1847	write(,) NF_STRERROR(ierr)
	1848	stop "getvarup"
	1849	endif
	1850	! write(,)'lecture T_adv_v ok'
	1851	!T adv in K/s
	1852	do l=1,ntime
	1853	do k=1,nlevel
	1854	T_adv_v(k,l)=T_adv_v(k,l)/3600.
	1855	enddo
	1856	enddo
	1857
	1858	#ifdef NC_DOUBLE
	1859	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
	1860	#else
	1861	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
	1862	#endif
	1863	if(ierr/=NF_NOERR) then
	1864	write(,) NF_STRERROR(ierr)
	1865	stop "getvarup"
	1866	endif
	1867	! write(,)'lecture q_adv_h ok'
	1868	!q adv in kg/kg/s
	1869	do l=1,ntime
	1870	do k=1,nlevel
	1871	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
	1872	enddo
	1873	enddo
	1874
	1875
	1876	#ifdef NC_DOUBLE
	1877	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
	1878	#else
	1879	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
	1880	#endif
	1881	if(ierr/=NF_NOERR) then
	1882	write(,) NF_STRERROR(ierr)
	1883	stop "getvarup"
	1884	endif
	1885	! write(,)'lecture q_adv_v ok'
	1886	!q adv in kg/kg/s
	1887	do l=1,ntime
	1888	do k=1,nlevel
	1889	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
	1890	enddo
	1891	enddo
	1892
	1893
	1894	#ifdef NC_DOUBLE
	1895	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
	1896	#else
	1897	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
	1898	#endif
	1899	if(ierr/=NF_NOERR) then
	1900	write(,) NF_STRERROR(ierr)
	1901	stop "getvarup"
	1902	endif
	1903
	1904	#ifdef NC_DOUBLE
	1905	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
	1906	#else
	1907	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
	1908	#endif
	1909	if(ierr/=NF_NOERR) then
	1910	write(,) NF_STRERROR(ierr)
	1911	stop "getvarup"
	1912	endif
	1913
	1914	#ifdef NC_DOUBLE
	1915	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
	1916	#else
	1917	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
	1918	#endif
	1919	if(ierr/=NF_NOERR) then
	1920	write(,) NF_STRERROR(ierr)
	1921	stop "getvarup"
	1922	endif
	1923
	1924	#ifdef NC_DOUBLE
	1925	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
	1926	#else
	1927	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
	1928	#endif
	1929	if(ierr/=NF_NOERR) then
	1930	write(,) NF_STRERROR(ierr)
	1931	stop "getvarup"
	1932	endif
	1933
	1934	#ifdef NC_DOUBLE
	1935	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
	1936	#else
	1937	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
	1938	#endif
	1939	if(ierr/=NF_NOERR) then
	1940	write(,) NF_STRERROR(ierr)
	1941	stop "getvarup"
	1942	endif
	1943
	1944	#ifdef NC_DOUBLE
	1945	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
	1946	#else
	1947	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
	1948	#endif
	1949	if(ierr/=NF_NOERR) then
	1950	write(,) NF_STRERROR(ierr)
	1951	stop "getvarup"
	1952	endif
	1953	! write(,)'lecture T_srf ok', T_srf
	1954
	1955	return
	1956	end subroutine read_twpice
	1957	!=====================================================================
	1958	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
	1959
	1960	! use netcdf
	1961
	1962	implicit none
	1963	#include "netcdf.inc"
	1964	integer nid,ttm,llm
	1965	real*8 time(ttm)
	1966	real*8 lev(llm)
	1967	integer ierr
	1968
	1969	integer i
	1970	integer timevar,levvar
	1971	integer timelen,levlen
	1972	integer timedimin,levdimin
	1973
	1974	! Control & lecture on dimensions
	1975	! ===============================
	1976	ierr=NF_INQ_DIMID(nid,"time",timedimin)
	1977	ierr=NF_INQ_VARID(nid,"time",timevar)
	1978	if (ierr.NE.NF_NOERR) then
	1979	write(,) 'ERROR: Field <time> is missing'
	1980	stop ""
	1981	endif
	1982	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
	1983
	1984	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
	1985	ierr=NF_INQ_VARID(nid,"lev",levvar)
	1986	if (ierr.NE.NF_NOERR) then
	1987	write(,) 'ERROR: Field <lev> is lacking'
	1988	stop ""
	1989	endif
	1990	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
	1991
	1992	if((timelen/=ttm).or.(levlen/=llm)) then
	1993	write(,) 'ERROR: Not the good lenght for axis'
	1994	write(,) 'longitude: ',timelen,ttm+1
	1995	write(,) 'latitude: ',levlen,llm
	1996	stop ""
	1997	endif
	1998
	1999	!#ifdef NC_DOUBLE
	2000	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
	2001	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
	2002	!#else
	2003	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
	2004	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
	2005	!#endif
	2006
	2007	return
	2008	end
	2009
	2010	!======================================================================
	2011	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play
	2012	: ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
	2013	: ,dth_dyn,dqh_dyn)
	2014	implicit none
	2015
	2016	c-------------------------------------------------------------------------
	2017	c Read RICO forcing data
	2018	c-------------------------------------------------------------------------
	2019	#include "dimensions.h"
	2020
	2021
	2022	integer nlev_rico
	2023	real ts_rico,ps_rico
	2024	real t_rico(llm),q_rico(llm)
	2025	real u_rico(llm),v_rico(llm)
	2026	real w_rico(llm)
	2027	real dth_dyn(llm)
	2028	real dqh_dyn(llm)
	2029
	2030
	2031	real play(llm),zlay(llm)
	2032
	2033
	2034	real prico(nlev_rico),zrico(nlev_rico)
	2035
	2036	character*80 fich_rico
	2037
	2038	integer k,l
	2039
	2040
	2041	print*,fich_rico
[1608]	2042	open(21,file=trim(fich_rico),form='formatted')
[1607]	2043	do k=1,llm
	2044	zlay(k)=0.
	2045	enddo
	2046
	2047	read(21,*) ps_rico,ts_rico
	2048	prico(1)=ps_rico
	2049	zrico(1)=0.0
	2050	do l=2,nlev_rico
	2051	read(21,*) k,prico(l),zrico(l)
	2052	enddo
	2053	close(21)
	2054
	2055	do k=1,llm
	2056	do l=1,80
	2057	if(prico(l)>play(k)) then
	2058	if(play(k)>prico(l+1)) then
	2059	zlay(k)=zrico(l)+(play(k)-prico(l)) *
	2060	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
	2061	else
	2062	zlay(k)=zrico(l)+(play(k)-prico(80))*
	2063	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
	2064	endif
	2065	endif
	2066	enddo
	2067	print*,k,zlay(k)
	2068	! U
	2069	if(0 < zlay(k) .and. zlay(k) < 4000) then
	2070	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
	2071	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
	2072	u_rico(k)= -1.9 + (30.0 + 1.9) /
	2073	: (12000 - 4000) * (zlay(k) - 4000)
	2074	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
	2075	u_rico(k)=30.0
	2076	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
	2077	u_rico(k)=30.0 - (30.0) /
	2078	: (20000 - 13000) * (zlay(k) - 13000)
	2079	else
	2080	u_rico(k)=0.0
	2081	endif
	2082
	2083	!Q_v
	2084	if(0 < zlay(k) .and. zlay(k) < 740) then
	2085	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
	2086	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
	2087	q_rico(k)=13.8 + (2.4 - 13.8) /
	2088	: (3260 - 740) * (zlay(k) - 740)
	2089	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
	2090	q_rico(k)=2.4 + (1.8 - 2.4) /
	2091	: (4000 - 3260) * (zlay(k) - 3260)
	2092	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
	2093	q_rico(k)=1.8 + (0 - 1.8) /
	2094	: (10000 - 4000) * (zlay(k) - 4000)
	2095	else
	2096	q_rico(k)=0.0
	2097	endif
	2098
	2099	!T
	2100	if(0 < zlay(k) .and. zlay(k) < 740) then
	2101	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
	2102	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
	2103	t_rico(k)=292.0 + (278.0 - 292.0) /
	2104	: (4000 - 740) * (zlay(k) - 740)
	2105	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
	2106	t_rico(k)=278.0 + (203.0 - 278.0) /
	2107	: (15000 - 4000) * (zlay(k) - 4000)
	2108	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
	2109	t_rico(k)=203.0 + (194.0 - 203.0) /
	2110	: (17500 - 15000)* (zlay(k) - 15000)
	2111	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
	2112	t_rico(k)=194.0 + (206.0 - 194.0) /
	2113	: (20000 - 17500)* (zlay(k) - 17500)
	2114	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
	2115	t_rico(k)=206.0 + (270.0 - 206.0) /
	2116	: (60000 - 20000)* (zlay(k) - 20000)
	2117	endif
	2118
	2119	! W
	2120	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
	2121	w_rico(k)=- (0.005/2260) * zlay(k)
	2122	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
	2123	w_rico(k)=- 0.005
	2124	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
	2125	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
	2126	else
	2127	w_rico(k)=0.0
	2128	endif
	2129
	2130	! dThrz+dTsw0+dTlw0
	2131	if(0 < zlay(k) .and. zlay(k) < 4000) then
	2132	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/
	2133	: (864004000) zlay(k)
	2134	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
	2135	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /
	2136	: (86400(5000 - 4000)) (zlay(k) - 4000)
	2137	else
	2138	dth_dyn(k)=0.0
	2139	endif
	2140	! dQhrz
	2141	if(0 < zlay(k) .and. zlay(k) < 3000) then
	2142	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/
	2143	: (864003000) (zlay(k))
	2144	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
	2145	dqh_dyn(k)=0.345 / 86400
	2146	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
	2147	dqh_dyn(k)=0.345 / 86400 +
	2148	+ (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
	2149	else
	2150	dqh_dyn(k)=0.0
	2151	endif
	2152
	2153	!? if(play(k)>6e4) then
	2154	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
	2155	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
	2156	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
	2157	!? *(6e4-play(k))/(6e4-3e4)
	2158	!? else
	2159	!? ratqs0(1,k)=ratqshaut
	2160	!? endif
	2161
	2162	enddo
	2163
	2164	do k=1,llm
	2165	q_rico(k)=q_rico(k)/1e3
	2166	dqh_dyn(k)=dqh_dyn(k)/1e3
	2167	v_rico(k)=-3.8
	2168	enddo
	2169
	2170	return
	2171	end
	2172
	2173	!=====================================================================
	2174	c-------------------------------------------------------------------------
	2175	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,
	2176	: sens,flat,adv_theta,rad_theta,adv_qt)
	2177	implicit none
	2178
	2179	c-------------------------------------------------------------------------
	2180	c Read ARM_CU case forcing data
	2181	c-------------------------------------------------------------------------
	2182
	2183	integer nlev_armcu,nt_armcu
	2184	real sens(nt_armcu),flat(nt_armcu)
	2185	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
	2186	character*80 fich_armcu
	2187
	2188	integer no,l,k,ip
	2189	real riy,rim,rid,rih,bid
	2190
	2191	integer iy,im,id,ih,in
	2192
[1608]	2193	open(21,file=trim(fich_armcu),form='formatted')
[1607]	2194	read(21,'(a)')
	2195	do ip = 1, nt_armcu
	2196	read(21,'(a)')
	2197	read(21,'(a)')
	2198	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),
	2199	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
	2200	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),
	2201	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
	2202	enddo
	2203	close(21)
	2204
	2205	223 format(5i3,5f8.3)
	2206	226 format(f7.1,1x,10f8.2)
	2207	227 format(f7.1,1x,1p,4e11.3)
	2208	230 format(6f9.3,4e11.3)
	2209
	2210	return
	2211	end
	2212
	2213	!=====================================================================
	2214	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof
	2215	: ,t_prof,q_prof,u_prof,v_prof,w_prof
	2216	: ,ht_prof,vt_prof,hq_prof,vq_prof
	2217	: ,t_mod,q_mod,u_mod,v_mod,w_mod
	2218	: ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
	2219
	2220	implicit none
	2221
	2222	#include "dimensions.h"
	2223
	2224	c-------------------------------------------------------------------------
	2225	c Vertical interpolation of TOGA-COARE forcing data onto model levels
	2226	c-------------------------------------------------------------------------
	2227
	2228	integer nlevmax
	2229	parameter (nlevmax=41)
	2230	integer nlev_toga,mxcalc
	2231	! real play(llm), plev_prof(nlevmax)
	2232	! real t_prof(nlevmax),q_prof(nlevmax)
	2233	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
	2234	! real ht_prof(nlevmax),vt_prof(nlevmax)
	2235	! real hq_prof(nlevmax),vq_prof(nlevmax)
	2236
	2237	real play(llm), plev_prof(nlev_toga)
	2238	real t_prof(nlev_toga),q_prof(nlev_toga)
	2239	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
	2240	real ht_prof(nlev_toga),vt_prof(nlev_toga)
	2241	real hq_prof(nlev_toga),vq_prof(nlev_toga)
	2242
	2243	real t_mod(llm),q_mod(llm)
	2244	real u_mod(llm),v_mod(llm), w_mod(llm)
	2245	real ht_mod(llm),vt_mod(llm)
	2246	real hq_mod(llm),vq_mod(llm)
	2247
	2248	integer l,k,k1,k2,kp
	2249	real aa,frac,frac1,frac2,fact
	2250
	2251	do l = 1, llm
	2252
	2253	if (play(l).ge.plev_prof(nlev_toga)) then
	2254
	2255	mxcalc=l
	2256	k1=0
	2257	k2=0
	2258
	2259	if (play(l).le.plev_prof(1)) then
	2260
	2261	do k = 1, nlev_toga-1
	2262	if (play(l).le.plev_prof(k)
	2263	: .and. play(l).gt.plev_prof(k+1)) then
	2264	k1=k
	2265	k2=k+1
	2266	endif
	2267	enddo
	2268
	2269	if (k1.eq.0 .or. k2.eq.0) then
	2270	write(,) 'PB! k1, k2 = ',k1,k2
	2271	write(,) 'l,play(l) = ',l,play(l)/100
	2272	do k = 1, nlev_toga-1
	2273	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
	2274	enddo
	2275	endif
	2276
	2277	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
	2278	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
	2279	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
	2280	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
	2281	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
	2282	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
	2283	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
	2284	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
	2285	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
	2286	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
	2287
	2288	else !play>plev_prof(1)
	2289
	2290	k1=1
	2291	k2=2
	2292	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
	2293	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
	2294	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
	2295	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
	2296	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
	2297	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
	2298	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
	2299	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
	2300	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
	2301	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
	2302	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
	2303
	2304	endif ! play.le.plev_prof(1)
	2305
	2306	else ! above max altitude of forcing file
	2307
	2308	cjyg
	2309	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
	2310	fact = max(fact,0.) !jyg
	2311	fact = exp(-fact) !jyg
	2312	t_mod(l)= t_prof(nlev_toga) !jyg
	2313	q_mod(l)= q_prof(nlev_toga)*fact !jyg
	2314	u_mod(l)= u_prof(nlev_toga)*fact !jyg
	2315	v_mod(l)= v_prof(nlev_toga)*fact !jyg
	2316	w_mod(l)= 0.0 !jyg
	2317	ht_mod(l)= ht_prof(nlev_toga) !jyg
	2318	vt_mod(l)= vt_prof(nlev_toga) !jyg
	2319	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
	2320	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
	2321
	2322	endif ! play
	2323
	2324	enddo ! l
	2325
	2326	! do l = 1,llm
	2327	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
	2328	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
	2329	! enddo
	2330
	2331	return
	2332	end
	2333
	2334	!======================================================================
	2335	SUBROUTINE interp_toga_time(day,day1,annee_ref
	2336	i ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga
	2337	i ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
	2338	i ,ht_toga,vt_toga,hq_toga,vq_toga
	2339	o ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
	2340	o ,ht_prof,vt_prof,hq_prof,vq_prof)
	2341	implicit none
	2342
	2343	!---------------------------------------------------------------------------------------
	2344	! Time interpolation of a 2D field to the timestep corresponding to day
	2345	!
	2346	! day: current julian day (e.g. 717538.2)
	2347	! day1: first day of the simulation
	2348	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
	2349	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
	2350	!---------------------------------------------------------------------------------------
	2351
	2352	#include "compar1d.h"
	2353
	2354	! inputs:
	2355	integer annee_ref
	2356	integer nt_toga,nlev_toga
	2357	integer year_ini_toga
	2358	real day, day1,day_ini_toga,dt_toga
	2359	real ts_toga(nt_toga)
	2360	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
	2361	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
	2362	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
	2363	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
	2364	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
	2365	! outputs:
	2366	real ts_prof
	2367	real plev_prof(nlev_toga),t_prof(nlev_toga)
	2368	real q_prof(nlev_toga),u_prof(nlev_toga)
	2369	real v_prof(nlev_toga),w_prof(nlev_toga)
	2370	real ht_prof(nlev_toga),vt_prof(nlev_toga)
	2371	real hq_prof(nlev_toga),vq_prof(nlev_toga)
	2372	! local:
	2373	integer it_toga1, it_toga2,k
	2374	real timeit,time_toga1,time_toga2,frac
	2375
	2376
	2377	if (forcing_type.eq.2) then
	2378	! Check that initial day of the simulation consistent with TOGA-COARE period:
	2379	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
	2380	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
	2381	print*,'Changer annee_ref dans run.def'
	2382	stop
	2383	endif
	2384	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
	2385	print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
	2386	print*,'Changer dayref dans run.def'
	2387	stop
	2388	endif
	2389	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
	2390	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
	2391	print*,'Changer dayref ou nday dans run.def'
	2392	stop
	2393	endif
	2394
	2395	else if (forcing_type.eq.4) then
	2396
	2397	! Check that initial day of the simulation consistent with TWP-ICE period:
	2398	if (annee_ref.ne.2006) then
	2399	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
	2400	print*,'Changer annee_ref dans run.def'
	2401	stop
	2402	endif
	2403	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
	2404	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
	2405	print*,'Changer dayref dans run.def'
	2406	stop
	2407	endif
	2408	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
	2409	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
	2410	print*,'Changer dayref ou nday dans run.def'
	2411	stop
	2412	endif
	2413
	2414	endif
	2415
	2416	! Determine timestep relative to the 1st day of TOGA-COARE:
	2417	! timeit=(day-day1)*86400.
	2418	! if (annee_ref.eq.1992) then
	2419	! timeit=(day-day_ini_toga)*86400.
	2420	! else
	2421	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
	2422	! endif
	2423	timeit=(day-day_ini_toga)*86400
	2424
	2425	! Determine the closest observation times:
	2426	it_toga1=INT(timeit/dt_toga)+1
	2427	it_toga2=it_toga1 + 1
	2428	time_toga1=(it_toga1-1)*dt_toga
	2429	time_toga2=(it_toga2-1)*dt_toga
	2430
	2431	if (it_toga1 .ge. nt_toga) then
	2432	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: '
	2433	: ,day,it_toga1,it_toga2,timeit/86400.
	2434	stop
	2435	endif
	2436
	2437	! time interpolation:
	2438	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
	2439	frac=max(frac,0.0)
	2440
	2441	ts_prof = ts_toga(it_toga2)
	2442	: -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
	2443
	2444	! print*,
	2445	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
	2446	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
	2447
	2448	do k=1,nlev_toga
	2449	plev_prof(k) = 100.*(plev_toga(k,it_toga2)
	2450	: -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
	2451	t_prof(k) = t_toga(k,it_toga2)
	2452	: -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
	2453	q_prof(k) = q_toga(k,it_toga2)
	2454	: -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
	2455	u_prof(k) = u_toga(k,it_toga2)
	2456	: -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
	2457	v_prof(k) = v_toga(k,it_toga2)
	2458	: -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
	2459	w_prof(k) = w_toga(k,it_toga2)
	2460	: -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
	2461	ht_prof(k) = ht_toga(k,it_toga2)
	2462	: -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
	2463	vt_prof(k) = vt_toga(k,it_toga2)
	2464	: -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
	2465	hq_prof(k) = hq_toga(k,it_toga2)
	2466	: -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
	2467	vq_prof(k) = vq_toga(k,it_toga2)
	2468	: -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
	2469	enddo
	2470
	2471	return
	2472	END
	2473
	2474	!======================================================================
	2475	SUBROUTINE interp_armcu_time(day,day1,annee_ref
	2476	i ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu
	2477	i ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu
	2478	i ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
	2479	implicit none
	2480
	2481	!---------------------------------------------------------------------------------------
	2482	! Time interpolation of a 2D field to the timestep corresponding to day
	2483	!
	2484	! day: current julian day (e.g. 717538.2)
	2485	! day1: first day of the simulation
	2486	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
	2487	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
	2488	! fs= sensible flux
	2489	! fl= latent flux
	2490	! at,rt,aqt= advective and radiative tendencies
	2491	!---------------------------------------------------------------------------------------
	2492
	2493	! inputs:
	2494	integer annee_ref
	2495	integer nt_armcu,nlev_armcu
	2496	integer year_ini_armcu
	2497	real day, day1,day_ini_armcu,dt_armcu
	2498	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
	2499	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
	2500	! outputs:
	2501	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
	2502	! local:
	2503	integer it_armcu1, it_armcu2,k
	2504	real timeit,time_armcu1,time_armcu2,frac
	2505
	2506	! Check that initial day of the simulation consistent with ARMCU period:
	2507	if (annee_ref.ne.1997 ) then
	2508	print*,'Pour ARMCU, annee_ref doit etre 1997 '
	2509	print*,'Changer annee_ref dans run.def'
	2510	stop
	2511	endif
	2512
	2513	timeit=(day-day_ini_armcu)*86400
	2514
	2515	! Determine the closest observation times:
	2516	it_armcu1=INT(timeit/dt_armcu)+1
	2517	it_armcu2=it_armcu1 + 1
	2518	time_armcu1=(it_armcu1-1)*dt_armcu
	2519	time_armcu2=(it_armcu2-1)*dt_armcu
	2520	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
	2521	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',
	2522	. it_armcu1,it_armcu2,time_armcu1,time_armcu2
	2523
	2524	if (it_armcu1 .ge. nt_armcu) then
	2525	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '
	2526	: ,day,it_armcu1,it_armcu2,timeit/86400.
	2527	stop
	2528	endif
	2529
	2530	! time interpolation:
	2531	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
	2532	frac=max(frac,0.0)
	2533
	2534	fs_prof = fs_armcu(it_armcu2)
	2535	: -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
	2536	fl_prof = fl_armcu(it_armcu2)
	2537	: -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
	2538	at_prof = at_armcu(it_armcu2)
	2539	: -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
	2540	rt_prof = rt_armcu(it_armcu2)
	2541	: -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
	2542	aqt_prof = aqt_armcu(it_armcu2)
	2543	: -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
	2544
	2545	print*,
	2546	:'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',
	2547	:day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,
	2548	:it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
	2549
	2550	return
	2551	END
	2552
	2553	!=====================================================================
	2554	subroutine readprofiles(nlev_max,kmax,height,
	2555	. thlprof,qtprof,uprof,
	2556	. vprof,e12prof,ugprof,vgprof,
	2557	. wfls,dqtdxls,dqtdyls,dqtdtls,
	2558	. thlpcar)
	2559	implicit none
	2560
	2561	integer nlev_max,kmax,kmax2
	2562	logical :: llesread = .true.
	2563
	2564	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),
	2565	. uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
	2566	. ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
	2567	. dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),
	2568	. thlpcar(nlev_max)
	2569
	2570	integer, parameter :: ilesfile=1
	2571	integer :: ierr,irad,imax,jtot,k
	2572	logical :: lmoist,lcoriol,ltimedep
	2573	real :: xsize,ysize
	2574	real :: ustin,wsvsurf,timerad
	2575	character(80) :: chmess
	2576
	2577	if(.not.(llesread)) return
	2578
	2579	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
	2580	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	2581	read (ilesfile,*) kmax
	2582	do k=1,kmax
	2583	read (ilesfile,*) height(k),thlprof(k),qtprof (k),
	2584	. uprof (k),vprof (k),e12prof(k)
	2585	enddo
	2586	close(ilesfile)
	2587
	2588	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
	2589	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
	2590	read (ilesfile,*) kmax2
	2591	if (kmax .ne. kmax2) then
	2592	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
	2593	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
	2594	stop 'lecture profiles'
	2595	endif
	2596	do k=1,kmax
	2597	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),
	2598	. dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
	2599	end do
	2600	close(ilesfile)
	2601
	2602	return
	2603	end
	2604
	2605
	2606	!======================================================================
	2607	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
	2608	. vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
	2609	!======================================================================
	2610	implicit none
	2611
	2612	integer nlev_max,kmax,kmax2
	2613	logical :: llesread = .true.
	2614
	2615	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
	2616	. thetaprof(nlev_max),rvprof(nlev_max),
	2617	. qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
	2618	. aprof(nlev_max+1),bprof(nlev_max+1)
	2619
	2620	integer, parameter :: ilesfile=1
	2621	integer, parameter :: ifile=2
	2622	integer :: ierr,irad,imax,jtot,k
	2623	logical :: lmoist,lcoriol,ltimedep
	2624	real :: xsize,ysize
	2625	real :: ustin,wsvsurf,timerad
	2626	character(80) :: chmess
	2627
	2628	if(.not.(llesread)) return
	2629
	2630	! Read profiles at full levels
	2631	IF(nlev_max.EQ.19) THEN
	2632	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
	2633	print *,'On ouvre prof.inp.19'
	2634	ELSE
	2635	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
	2636	print *,'On ouvre prof.inp.40'
	2637	ENDIF
	2638	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
	2639	read (ilesfile,*) kmax
	2640	do k=1,kmax
	2641	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k),
	2642	. thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
	2643	enddo
	2644	close(ilesfile)
	2645
	2646	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
	2647	IF(nlev_max.EQ.19) THEN
	2648	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
	2649	print *,'On ouvre proh.inp.19'
	2650	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
	2651	ELSE
	2652	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
	2653	print *,'On ouvre proh.inp.40'
	2654	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
	2655	ENDIF
	2656	read (ifile,*) kmax
	2657	do k=1,kmax
	2658	read (ifile,*) jtot,aprof(k),bprof(k)
	2659	enddo
	2660	close(ifile)
	2661
	2662	return
	2663	end
	2664	!===============================================================
	2665	function ismin(n,sx,incx)
	2666
	2667	implicit none
	2668	integer n,i,incx,ismin,ix
	2669	real sx((n-1)*incx+1),sxmin
	2670
	2671	ix=1
	2672	ismin=1
	2673	sxmin=sx(1)
	2674	do i=1,n-1
	2675	ix=ix+incx
	2676	if(sx(ix).lt.sxmin) then
	2677	sxmin=sx(ix)
	2678	ismin=i+1
	2679	endif
	2680	enddo
	2681
	2682	return
	2683	end
	2684
	2685	!===============================================================
	2686	function ismax(n,sx,incx)
	2687
	2688	implicit none
	2689	integer n,i,incx,ismax,ix
	2690	real sx((n-1)*incx+1),sxmax
	2691
	2692	ix=1
	2693	ismax=1
	2694	sxmax=sx(1)
	2695	do i=1,n-1
	2696	ix=ix+incx
	2697	if(sx(ix).gt.sxmax) then
	2698	sxmax=sx(ix)
	2699	ismax=i+1
	2700	endif
	2701	enddo
	2702
	2703	return
	2704	end
	2705

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

Download in other formats:

Original Format