Context Navigation

1DUTILS.h_with_writelim_old @ 1794

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

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

Download in other formats:

Original Format