Context Navigation

guide_mod.F90 @ 1721

Last change on this file since 1721 was 1720, checked in by Ehouarn Millour, 12 years ago
Some cleanup around pres2lev: the three dynamical cores now use the module pres2lev_mod.F90 EM
Property svn:executable set to ``*
File size: 52.0 KB

Rev	Line
[1170]	1	!
[1186]	2	! $Id$
[1170]	3	!
	4	MODULE guide_mod
	5
	6	!=======================================================================
	7	! Auteur: F.Hourdin
	8	! F. Codron 01/09
	9	!=======================================================================
	10
	11	USE getparam
	12	USE Write_Field
[1188]	13	use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
[1720]	14	use pres2lev_mod
[1170]	15
	16	IMPLICIT NONE
	17
	18	! ---------------------------------------------
	19	! Declarations des cles logiques et parametres
	20	! ---------------------------------------------
	21	INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
	22	INTEGER, PRIVATE, SAVE :: nlevnc
	23	LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P
	24	LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta
	25	LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
	26	LOGICAL, PRIVATE, SAVE :: guide_modele,invert_p,invert_y,ini_anal
	27	LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav
	28
	29	REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
	30	REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
	31	REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
	32	REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
	33	REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
	34
	35	REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
	36	REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
	37	REAL, PRIVATE, SAVE :: tau_lon,tau_lat
	38
	39	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
	40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q
	41	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
	42
	43	! ---------------------------------------------
	44	! Variables de guidage
	45	! ---------------------------------------------
	46	! Variables des fichiers de guidage
	47	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
	48	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
	49	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
	50	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
	51	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2
	52	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc
	53	! Variables aux dimensions du modele
	54	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2
	55	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2
	56	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2
	57	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2
	58	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2
	59
	60	CONTAINS
	61	!=======================================================================
	62
	63	SUBROUTINE guide_init
	64
[1403]	65	USE control_mod
	66
[1170]	67	IMPLICIT NONE
	68
	69	INCLUDE "dimensions.h"
	70	INCLUDE "paramet.h"
	71	INCLUDE "netcdf.inc"
	72
	73	INTEGER :: error,ncidpl,rid,rcod
	74	CHARACTER (len = 80) :: abort_message
	75	CHARACTER (len = 20) :: modname = 'guide_init'
	76
	77	! ---------------------------------------------
	78	! Lecture des parametres:
	79	! ---------------------------------------------
	80	! Variables guidees
	81	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
	82	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
	83	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
	84	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
	85	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
	86	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
	87	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
	88
	89	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
	90	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
	91
	92	! Constantes de rappel. Unite : fraction de jour
	93	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
	94	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
	95	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
	96	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
	97	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
	98	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
	99	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
	100	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
	101	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
	102	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
	103	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
	104	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
	105
	106	! Sauvegarde du for�age
	107	CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
	108	CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
	109	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
	110	IF (iguide_sav.GT.0) THEN
	111	iguide_sav=day_step/iguide_sav
	112	ELSE
	113	iguide_sav=day_step*iguide_sav
	114	ENDIF
	115
	116	! Guidage regional seulement (sinon constant ou suivant le zoom)
	117	CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
	118	CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
	119	CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
	120	CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
	121	CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
	122	CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
	123	CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
	124
	125	! Parametres pour lecture des fichiers
	126	CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
	127	CALL getpar('iguide_int',4,iguide_int,'freq. lecture guidage')
	128	IF (iguide_int.GT.0) THEN
	129	iguide_int=day_step/iguide_int
	130	ELSE
	131	iguide_int=day_step*iguide_int
	132	ENDIF
	133	CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele')
	134	CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
	135	CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
	136	CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
	137	CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
	138
	139	! ---------------------------------------------
	140	! Determination du nombre de niveaux verticaux
	141	! des fichiers guidage
	142	! ---------------------------------------------
	143	ncidpl=-99
	144	if (guide_modele) then
[1188]	145	if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
[1170]	146	else
	147	if (guide_u) then
[1188]	148	if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
[1170]	149	elseif (guide_v) then
[1188]	150	if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
[1170]	151	elseif (guide_T) then
[1188]	152	if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
[1170]	153	elseif (guide_Q) then
[1188]	154	if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
[1170]	155	endif
	156	endif
	157	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
	158	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
	159	IF (error.NE.NF_NOERR) THEN
	160	print *,'Guide: probleme lecture niveaux pression'
	161	CALL abort_gcm(modname,abort_message,1)
	162	ENDIF
	163	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
	164	print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
[1188]	165	rcod = nf90_close(ncidpl)
[1170]	166
	167	! ---------------------------------------------
	168	! Allocation des variables
	169	! ---------------------------------------------
	170	abort_message='pb in allocation guide'
	171
	172	ALLOCATE(apnc(nlevnc), stat = error)
	173	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	174	ALLOCATE(bpnc(nlevnc), stat = error)
	175	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	176	apnc=0.;bpnc=0.
	177
	178	ALLOCATE(alpha_pcor(llm), stat = error)
	179	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	180	ALLOCATE(alpha_u(ip1jmp1), stat = error)
	181	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	182	ALLOCATE(alpha_v(ip1jm), stat = error)
	183	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	184	ALLOCATE(alpha_T(ip1jmp1), stat = error)
	185	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	186	ALLOCATE(alpha_Q(ip1jmp1), stat = error)
	187	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	188	ALLOCATE(alpha_P(ip1jmp1), stat = error)
	189	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	190	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
	191
	192	IF (guide_u) THEN
	193	ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
	194	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	195	ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
	196	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	197	ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
	198	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	199	ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
	200	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	201	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
	202	ENDIF
	203
	204	IF (guide_T) THEN
	205	ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
	206	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	207	ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
	208	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	209	ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
	210	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	211	ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
	212	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	213	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
	214	ENDIF
	215
	216	IF (guide_Q) THEN
	217	ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
	218	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	219	ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
	220	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	221	ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
	222	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	223	ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
	224	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	225	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
	226	ENDIF
	227
	228	IF (guide_v) THEN
	229	ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
	230	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	231	ALLOCATE(vgui1(ip1jm,llm), stat = error)
	232	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	233	ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
	234	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	235	ALLOCATE(vgui2(ip1jm,llm), stat = error)
	236	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	237	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
	238	ENDIF
	239
	240	IF (guide_P.OR.guide_modele) THEN
	241	ALLOCATE(psnat1(iip1,jjp1), stat = error)
	242	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	243	ALLOCATE(psnat2(iip1,jjp1), stat = error)
	244	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	245	psnat1=0.;psnat2=0.;
	246	ENDIF
	247	IF (guide_P) THEN
	248	ALLOCATE(psgui2(ip1jmp1), stat = error)
	249	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	250	ALLOCATE(psgui1(ip1jmp1), stat = error)
	251	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	252	psgui1=0.;psgui2=0.
	253	ENDIF
	254
	255	! ---------------------------------------------
	256	! Lecture du premier etat de guidage.
	257	! ---------------------------------------------
	258	IF (guide_2D) THEN
	259	CALL guide_read2D(1)
	260	ELSE
	261	CALL guide_read(1)
	262	ENDIF
	263	IF (guide_v) vnat1=vnat2
	264	IF (guide_u) unat1=unat2
	265	IF (guide_T) tnat1=tnat2
	266	IF (guide_Q) qnat1=qnat2
	267	IF (guide_P.OR.guide_modele) psnat1=psnat2
	268
	269	END SUBROUTINE guide_init
	270
	271	!=======================================================================
	272	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
[1403]	273
	274	USE control_mod
[1170]	275
	276	IMPLICIT NONE
	277
	278	INCLUDE "dimensions.h"
	279	INCLUDE "paramet.h"
	280	INCLUDE "comconst.h"
	281	INCLUDE "comvert.h"
	282
	283	! Variables entree
	284	INTEGER, INTENT(IN) :: itau !pas de temps
	285	REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
	286	REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
	287	REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
	288
	289	! Variables locales
	290	LOGICAL, SAVE :: first=.TRUE.
	291	LOGICAL :: f_out ! sortie guidage
	292	REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
	293	REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P
	294	! Compteurs temps:
	295	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
	296	REAL :: ditau, dday_step
	297	REAL :: tau,reste ! position entre 2 etats de guidage
	298	REAL, SAVE :: factt ! pas de temps en fraction de jour
	299
	300	INTEGER :: l
	301
	302	!-----------------------------------------------------------------------
	303	! Initialisations au premier passage
	304	!-----------------------------------------------------------------------
	305	IF (first) THEN
	306	first=.FALSE.
	307	CALL guide_init
	308	itau_test=1001
	309	step_rea=1
	310	count_no_rea=0
	311	! Calcul des constantes de rappel
	312	factt=dtvr*iperiod/daysec
	313	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
	314	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
	315	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
	316	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
	317	call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
	318	! correction de rappel dans couche limite
	319	if (guide_BL) then
	320	alpha_pcor(:)=1.
	321	else
	322	do l=1,llm
	323	alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
	324	enddo
	325	endif
	326	! ini_anal: etat initial egal au guidage
	327	IF (ini_anal) THEN
	328	CALL guide_interp(ps,teta)
	329	IF (guide_u) ucov=ugui2
	330	IF (guide_v) vcov=ugui2
	331	IF (guide_T) teta=tgui2
	332	IF (guide_Q) q=qgui2
	333	IF (guide_P) THEN
	334	ps=psgui2
	335	CALL pression(ip1jmp1,ap,bp,ps,p)
	336	CALL massdair(p,masse)
	337	ENDIF
	338	RETURN
	339	ENDIF
	340	! Verification structure guidage
	341	IF (guide_u) THEN
	342	CALL writefield('unat',unat1)
	343	CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
	344	ENDIF
	345	IF (guide_T) THEN
	346	CALL writefield('tnat',tnat1)
	347	CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
	348	ENDIF
	349
	350	ENDIF !first
	351
	352	!-----------------------------------------------------------------------
	353	! Lecture des fichiers de guidage ?
	354	!-----------------------------------------------------------------------
	355	IF (iguide_read.NE.0) THEN
	356	ditau=real(itau)
	357	dday_step=real(day_step)
	358	IF (iguide_read.LT.0) THEN
[1403]	359	tau=ditau/dday_step/REAL(iguide_read)
[1170]	360	ELSE
[1403]	361	tau=REAL(iguide_read)*ditau/dday_step
[1170]	362	ENDIF
	363	reste=tau-AINT(tau)
	364	IF (reste.EQ.0.) THEN
	365	IF (itau_test.EQ.itau) THEN
	366	write(,)'deuxieme passage de advreel a itau=',itau
	367	stop
	368	ELSE
	369	IF (guide_v) vnat1=vnat2
	370	IF (guide_u) unat1=unat2
	371	IF (guide_T) tnat1=tnat2
	372	IF (guide_Q) qnat1=qnat2
	373	IF (guide_P.OR.guide_modele) psnat1=psnat2
	374	step_rea=step_rea+1
	375	itau_test=itau
	376	print*,'Lecture fichiers guidage, pas ',step_rea, &
	377	'apres ',count_no_rea,' non lectures'
	378	IF (guide_2D) THEN
	379	CALL guide_read2D(step_rea)
	380	ELSE
	381	CALL guide_read(step_rea)
	382	ENDIF
	383	count_no_rea=0
	384	ENDIF
	385	ELSE
	386	count_no_rea=count_no_rea+1
	387
	388	ENDIF
	389	ENDIF !iguide_read=0
	390
	391	!-----------------------------------------------------------------------
	392	! Interpolation et conversion des champs de guidage
	393	!-----------------------------------------------------------------------
	394	IF (MOD(itau,iguide_int).EQ.0) THEN
	395	CALL guide_interp(ps,teta)
	396	ENDIF
	397	! Repartition entre 2 etats de guidage
	398	IF (iguide_read.NE.0) THEN
	399	tau=reste
	400	ELSE
	401	tau=1.
	402	ENDIF
	403
	404	!-----------------------------------------------------------------------
	405	! Ajout des champs de guidage
	406	!-----------------------------------------------------------------------
	407	! Sauvegarde du guidage?
	408	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
	409	IF (f_out) CALL guide_out("S",jjp1,1,ps)
	410
	411	if (guide_u) then
	412	if (guide_add) then
	413	f_add=(1.-tau)ugui1+tauugui2
	414	else
	415	f_add=(1.-tau)ugui1+tauugui2-ucov
	416	endif
	417	if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
	418	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
	419	IF (f_out) CALL guide_out("U",jjp1,llm,f_add/factt)
	420	ucov=ucov+f_add
	421	endif
	422
	423	if (guide_T) then
	424	if (guide_add) then
	425	f_add=(1.-tau)tgui1+tautgui2
	426	else
	427	f_add=(1.-tau)tgui1+tautgui2-teta
	428	endif
	429	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
	430	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
	431	IF (f_out) CALL guide_out("T",jjp1,llm,f_add/factt)
	432	teta=teta+f_add
	433	endif
	434
	435	if (guide_P) then
	436	if (guide_add) then
	437	f_add(1:ip1jmp1,1)=(1.-tau)psgui1+taupsgui2
	438	else
	439	f_add(1:ip1jmp1,1)=(1.-tau)psgui1+taupsgui2-ps
	440	endif
	441	if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
	442	CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
	443	IF (f_out) CALL guide_out("P",jjp1,1,f_add(1:ip1jmp1,1)/factt)
	444	ps=ps+f_add(1:ip1jmp1,1)
	445	CALL pression(ip1jmp1,ap,bp,ps,p)
	446	CALL massdair(p,masse)
	447	endif
	448
	449	if (guide_Q) then
	450	if (guide_add) then
	451	f_add=(1.-tau)qgui1+tauqgui2
	452	else
	453	f_add=(1.-tau)qgui1+tauqgui2-q
	454	endif
	455	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
	456	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
	457	IF (f_out) CALL guide_out("Q",jjp1,llm,f_add/factt)
	458	q=q+f_add
	459	endif
	460
	461	if (guide_v) then
	462	if (guide_add) then
	463	f_add(1:ip1jm,:)=(1.-tau)vgui1+tauvgui2
	464	else
	465	f_add(1:ip1jm,:)=(1.-tau)vgui1+tauvgui2-vcov
	466	endif
	467	if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
	468	CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
	469	IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:)/factt)
	470	vcov=vcov+f_add(1:ip1jm,:)
	471	endif
	472	END SUBROUTINE guide_main
	473
	474	!=======================================================================
	475	SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
	476	! field1=afield1+alphafield2
	477
	478	IMPLICIT NONE
	479
	480	! input variables
	481	INTEGER, INTENT(IN) :: hsize
	482	INTEGER, INTENT(IN) :: vsize
	483	REAL, DIMENSION(hsize), INTENT(IN) :: alpha
	484	REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
	485
	486	! Local variables
	487	INTEGER :: l
	488
	489	do l=1,vsize
	490	field(:,l)=alphafield(:,l)alpha_pcor(l)
	491	enddo
	492
	493	END SUBROUTINE guide_addfield
	494
	495	!=======================================================================
	496	SUBROUTINE guide_zonave(typ,hsize,vsize,field)
	497
	498	IMPLICIT NONE
	499
	500	INCLUDE "dimensions.h"
	501	INCLUDE "paramet.h"
	502	INCLUDE "comgeom.h"
	503	INCLUDE "comconst.h"
	504
	505	! input/output variables
	506	INTEGER, INTENT(IN) :: typ
	507	INTEGER, INTENT(IN) :: vsize
	508	INTEGER, INTENT(IN) :: hsize
	509	REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
	510
	511	! Local variables
	512	LOGICAL, SAVE :: first=.TRUE.
	513	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
	514	INTEGER :: i,j,l,ij
	515	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
	516	REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
	517
	518	IF (first) THEN
	519	first=.FALSE.
	520	!Compute domain for averaging
	521	lond=rlonu*180./pi
	522	imin(1)=1;imax(1)=iip1;
	523	imin(2)=1;imax(2)=iip1;
	524	IF (guide_reg) THEN
	525	DO i=1,iim
	526	IF (lond(i).LT.lon_min_g) imin(1)=i
	527	IF (lond(i).LE.lon_max_g) imax(1)=i
	528	ENDDO
	529	lond=rlonv*180./pi
	530	DO i=1,iim
	531	IF (lond(i).LT.lon_min_g) imin(2)=i
	532	IF (lond(i).LE.lon_max_g) imax(2)=i
	533	ENDDO
	534	ENDIF
	535	ENDIF
	536
	537	fieldm=0.
	538	DO l=1,vsize
	539	! Compute zonal average
	540	DO j=1,hsize
	541	DO i=imin(typ),imax(typ)
	542	ij=(j-1)*iip1+i
	543	fieldm(j,l)=fieldm(j,l)+field(ij,l)
	544	ENDDO
	545	ENDDO
[1403]	546	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
[1170]	547	! Compute forcing
	548	DO j=1,hsize
	549	DO i=1,iip1
	550	ij=(j-1)*iip1+i
	551	field(ij,l)=fieldm(j,l)
	552	ENDDO
	553	ENDDO
	554	ENDDO
	555
	556	END SUBROUTINE guide_zonave
	557
	558	!=======================================================================
	559	SUBROUTINE guide_interp(psi,teta)
	560
	561	IMPLICIT NONE
	562
	563	include "dimensions.h"
	564	include "paramet.h"
	565	include "comvert.h"
	566	include "comgeom2.h"
	567	include "comconst.h"
	568
	569	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
	570	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
	571
	572	LOGICAL, SAVE :: first=.TRUE.
	573	! Variables pour niveaux pression:
	574	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
	575	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
	576	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
	577	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
	578	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
	579	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
	580	REAL, DIMENSION (iip1,jjm,llm) :: pbary
	581	! Variables pour fonction Exner (P milieu couche)
	582	REAL, DIMENSION (iip1,jjp1,llm) :: pk, pkf
	583	REAL, DIMENSION (iip1,jjp1,llm) :: alpha, beta
	584	REAL, DIMENSION (iip1,jjp1) :: pks
	585	REAL :: prefkap,unskap
	586	! Pression de vapeur saturante
	587	REAL, DIMENSION (ip1jmp1,llm) :: qsat
	588	!Variables intermediaires interpolation
	589	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
	590	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
	591
	592	INTEGER :: i,j,l,ij
	593
	594	print *,'Guide: conversion variables guidage'
	595	! -----------------------------------------------------------------
	596	! Calcul des niveaux de pression champs guidage
	597	! -----------------------------------------------------------------
	598	if (guide_modele) then
	599	do i=1,iip1
	600	do j=1,jjp1
	601	do l=1,nlevnc
	602	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
	603	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
	604	enddo
	605	enddo
	606	enddo
	607	else
	608	do i=1,iip1
	609	do j=1,jjp1
	610	do l=1,nlevnc
	611	plnc2(i,j,l)=apnc(l)
	612	plnc1(i,j,l)=apnc(l)
	613	enddo
	614	enddo
	615	enddo
	616
	617	endif
	618	if (first) then
	619	first=.FALSE.
	620	print*,'Guide: verification ordre niveaux verticaux'
	621	print*,'LMDZ :'
	622	do l=1,llm
	623	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
	624	+psi(1,jjp1)*(bp(l)+bp(l+1))/2.
	625	enddo
	626	print*,'Fichiers guidage'
	627	do l=1,nlevnc
	628	print*,'PL(',l,')=',plnc2(1,1,l)
	629	enddo
	630	print *,'inversion de l''ordre: invert_p=',invert_p
	631	if (guide_u) then
	632	do l=1,nlevnc
	633	print*,'U(',l,')=',unat2(1,1,l)
	634	enddo
	635	endif
	636	if (guide_T) then
	637	do l=1,nlevnc
	638	print*,'T(',l,')=',tnat2(1,1,l)
	639	enddo
	640	endif
	641	endif
	642
	643	! -----------------------------------------------------------------
	644	! Calcul niveaux pression modele
	645	! -----------------------------------------------------------------
	646	CALL pression( ip1jmp1, ap, bp, psi, p )
[1625]	647	if (pressure_exner) then
[1520]	648	CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
[1625]	649	else
[1520]	650	CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
	651	endif
[1170]	652	! .... Calcul de pls , pression au milieu des couches ,en Pascals
	653	unskap=1./kappa
	654	prefkap = preff ** kappa
	655	DO l = 1, llm
	656	DO j=1,jjp1
	657	DO i =1, iip1
	658	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
	659	ENDDO
	660	ENDDO
	661	ENDDO
	662
	663	! calcul des pressions pour les grilles u et v
	664	do l=1,llm
	665	do j=1,jjp1
	666	do i=1,iip1
	667	pext(i,j,l)=pls(i,j,l)*aire(i,j)
	668	enddo
	669	enddo
	670	enddo
	671	call massbar(pext, pbarx, pbary )
	672	do l=1,llm
	673	do j=1,jjp1
	674	do i=1,iip1
	675	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
	676	plsnc(i,j,l)=pls(i,j,l)
	677	enddo
	678	enddo
	679	enddo
	680	do l=1,llm
	681	do j=1,jjm
	682	do i=1,iip1
	683	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
	684	enddo
	685	enddo
	686	enddo
	687
	688	! -----------------------------------------------------------------
	689	! Interpolation champs guidage sur niveaux modele (+inversion N/S)
	690	! Conversion en variables gcm (ucov, vcov...)
	691	! -----------------------------------------------------------------
	692	if (guide_P) then
	693	do j=1,jjp1
	694	do i=1,iim
	695	ij=(j-1)*iip1+i
	696	psgui1(ij)=psnat1(i,j)
	697	psgui2(ij)=psnat2(i,j)
	698	enddo
	699	psgui1(iip1*j)=psnat1(1,j)
	700	psgui2(iip1*j)=psnat2(1,j)
	701	enddo
	702	endif
	703
	704	IF (guide_u) THEN
	705	CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
	706	CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
	707	do l=1,llm
	708	do j=1,jjp1
	709	do i=1,iim
	710	ij=(j-1)*iip1+i
	711	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
	712	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
	713	enddo
	714	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
	715	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
	716	enddo
	717	do i=1,iip1
	718	ugui1(i,l)=0.
	719	ugui1(ip1jm+i,l)=0.
	720	ugui2(i,l)=0.
	721	ugui2(ip1jm+i,l)=0.
	722	enddo
	723	enddo
	724	ENDIF
	725
	726	IF (guide_T) THEN
	727	CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
	728	CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
	729	do l=1,llm
	730	do j=1,jjp1
	731	IF (guide_teta) THEN
	732	do i=1,iim
	733	ij=(j-1)*iip1+i
	734	tgui1(ij,l)=zu1(i,j,l)
	735	tgui2(ij,l)=zu2(i,j,l)
	736	enddo
	737	ELSE
	738	do i=1,iim
	739	ij=(j-1)*iip1+i
	740	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
	741	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
	742	enddo
	743	ENDIF
	744	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
	745	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
	746	enddo
	747	do i=1,iip1
	748	tgui1(i,l)=tgui1(1,l)
	749	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
	750	tgui2(i,l)=tgui2(1,l)
	751	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
	752	enddo
	753	enddo
	754	ENDIF
	755
	756	IF (guide_v) THEN
	757
	758	CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
	759	CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
	760
	761	do l=1,llm
	762	do j=1,jjm
	763	do i=1,iim
	764	ij=(j-1)*iip1+i
	765	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
	766	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
	767	enddo
	768	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
	769	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
	770	enddo
	771	enddo
	772	ENDIF
	773
	774	IF (guide_Q) THEN
	775	! On suppose qu'on a la bonne variable dans le fichier de guidage:
	776	! Hum.Rel si guide_hr, Hum.Spec. sinon.
	777	CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
	778	CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
	779	do l=1,llm
	780	do j=1,jjp1
	781	do i=1,iim
	782	ij=(j-1)*iip1+i
	783	qgui1(ij,l)=zu1(i,j,l)
	784	qgui2(ij,l)=zu2(i,j,l)
	785	enddo
	786	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
	787	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
	788	enddo
	789	do i=1,iip1
	790	qgui1(i,l)=qgui1(1,l)
	791	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
	792	qgui2(i,l)=qgui2(1,l)
	793	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
	794	enddo
	795	enddo
	796	IF (guide_hr) THEN
	797	CALL q_sat(iip1jjp1llm,teta*pk/cpp,plsnc,qsat)
	798	qgui1=qgui1qsat0.01 !hum. rel. en %
	799	qgui2=qgui2qsat0.01
	800	ENDIF
	801	ENDIF
	802
	803	END SUBROUTINE guide_interp
	804
	805	!=======================================================================
	806	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
	807
	808	! Calcul des constantes de rappel alpha (=1/tau)
	809
	810	implicit none
	811
	812	include "dimensions.h"
	813	include "paramet.h"
	814	include "comconst.h"
	815	include "comgeom2.h"
	816	include "serre.h"
	817
	818	! input arguments :
	819	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
	820	INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
	821	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
	822	REAL, INTENT(IN) :: taumin,taumax
	823	! output arguments:
	824	REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
	825
	826	! local variables:
	827	LOGICAL, SAVE :: first=.TRUE.
	828	REAL, SAVE :: gamma,dxdy_min,dxdy_max
	829	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
	830	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
	831	REAL, DIMENSION (iip1,jjm) :: dxdyv
	832	real dxdy_
	833	real zlat,zlon
	834	real alphamin,alphamax,xi
	835	integer i,j,ilon,ilat
	836
	837
	838	alphamin=factt/taumax
	839	alphamax=factt/taumin
	840	IF (guide_reg.OR.guide_add) THEN
	841	alpha=alphamax
	842	!-----------------------------------------------------------------------
	843	! guide_reg: alpha=alpha_min dans region, 0. sinon.
	844	!-----------------------------------------------------------------------
	845	IF (guide_reg) THEN
	846	do j=1,pjm
	847	do i=1,pim
	848	if (typ.eq.2) then
	849	zlat=rlatu(j)*180./pi
	850	zlon=rlonu(i)*180./pi
	851	elseif (typ.eq.1) then
	852	zlat=rlatu(j)*180./pi
	853	zlon=rlonv(i)*180./pi
	854	elseif (typ.eq.3) then
	855	zlat=rlatv(j)*180./pi
	856	zlon=rlonv(i)*180./pi
	857	endif
	858	alpha(i,j)=alphamax/16.* &
	859	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
	860	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
	861	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
	862	(1.+tanh((lon_max_g-zlon)/tau_lon))
	863	enddo
	864	enddo
	865	ENDIF
	866	ELSE
	867	!-----------------------------------------------------------------------
	868	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
	869	!-----------------------------------------------------------------------
	870	!Calcul de l'aire des mailles
	871	do j=2,jjm
	872	do i=2,iip1
	873	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
	874	enddo
	875	zdx(1,j)=zdx(iip1,j)
	876	enddo
	877	do j=2,jjm
	878	do i=1,iip1
	879	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
	880	enddo
	881	enddo
	882	do i=1,iip1
	883	zdx(i,1)=zdx(i,2)
	884	zdx(i,jjp1)=zdx(i,jjm)
	885	zdy(i,1)=zdy(i,2)
	886	zdy(i,jjp1)=zdy(i,jjm)
	887	enddo
	888	do j=1,jjp1
	889	do i=1,iip1
	890	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
	891	enddo
	892	enddo
	893	IF (typ.EQ.2) THEN
	894	do j=1,jjp1
	895	do i=1,iim
	896	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
	897	enddo
	898	dxdyu(iip1,j)=dxdyu(1,j)
	899	enddo
	900	ENDIF
	901	IF (typ.EQ.3) THEN
	902	do j=1,jjm
	903	do i=1,iip1
	904	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
	905	enddo
	906	enddo
	907	ENDIF
	908	! Premier appel: calcul des aires min et max et de gamma.
	909	IF (first) THEN
	910	first=.FALSE.
	911	! coordonnees du centre du zoom
	912	CALL coordij(clon,clat,ilon,ilat)
	913	! aire de la maille au centre du zoom
	914	dxdy_min=dxdys(ilon,ilat)
	915	! dxdy maximale de la maille
	916	dxdy_max=0.
	917	do j=1,jjp1
	918	do i=1,iip1
	919	dxdy_max=max(dxdy_max,dxdys(i,j))
	920	enddo
	921	enddo
	922	! Calcul de gamma
	923	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
	924	print*,'ATTENTION modele peu zoome'
	925	print*,'ATTENTION on prend une constante de guidage cste'
	926	gamma=0.
	927	else
	928	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
	929	print*,'gamma=',gamma
	930	if (gamma.lt.1.e-5) then
	931	print*,'gamma =',gamma,'<1e-5'
	932	stop
	933	endif
	934	gamma=log(0.5)/log(gamma)
	935	if (gamma4) then
	936	gamma=min(gamma,4.)
	937	endif
	938	print*,'gamma=',gamma
	939	endif
	940	ENDIF !first
	941
	942	do j=1,pjm
	943	do i=1,pim
	944	if (typ.eq.1) then
	945	dxdy_=dxdys(i,j)
	946	zlat=rlatu(j)*180./pi
	947	elseif (typ.eq.2) then
	948	dxdy_=dxdyu(i,j)
	949	zlat=rlatu(j)*180./pi
	950	elseif (typ.eq.3) then
	951	dxdy_=dxdyv(i,j)
	952	zlat=rlatv(j)*180./pi
	953	endif
	954	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
	955	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
	956	alpha(i,j)=alphamin
	957	else
	958	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
	959	xi=min(xi,1.)
	960	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
	961	alpha(i,j)=xialphamin+(1.-xi)alphamax
	962	else
	963	alpha(i,j)=0.
	964	endif
	965	endif
	966	enddo
	967	enddo
	968	ENDIF ! guide_reg
	969
	970	END SUBROUTINE tau2alpha
	971
	972	!=======================================================================
	973	SUBROUTINE guide_read(timestep)
	974
	975	IMPLICIT NONE
	976
	977	#include "netcdf.inc"
	978	#include "dimensions.h"
	979	#include "paramet.h"
	980
	981	INTEGER, INTENT(IN) :: timestep
	982
	983	LOGICAL, SAVE :: first=.TRUE.
	984	! Identification fichiers et variables NetCDF:
	985	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
	986	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
	987	INTEGER :: ncidpl,varidpl,varidap,varidbp
	988	! Variables auxiliaires NetCDF:
	989	INTEGER, DIMENSION(4) :: start,count
	990	INTEGER :: status,rcode
	991
	992	! -----------------------------------------------------------------
	993	! Premier appel: initialisation de la lecture des fichiers
	994	! -----------------------------------------------------------------
	995	if (first) then
	996	ncidpl=-99
	997	print*,'Guide: ouverture des fichiers guidage '
	998	! Niveaux de pression si non constants
	999	if (guide_modele) then
	1000	print *,'Lecture du guidage sur niveaux mod�le'
[1188]	1001	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
	1002	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
	1003	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[1170]	1004	print*,'ncidpl,varidap',ncidpl,varidap
	1005	endif
	1006	! Vent zonal
	1007	if (guide_u) then
[1188]	1008	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
	1009	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[1170]	1010	print*,'ncidu,varidu',ncidu,varidu
	1011	if (ncidpl.eq.-99) ncidpl=ncidu
	1012	endif
	1013	! Vent meridien
	1014	if (guide_v) then
[1188]	1015	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
	1016	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[1170]	1017	print*,'ncidv,varidv',ncidv,varidv
	1018	if (ncidpl.eq.-99) ncidpl=ncidv
	1019	endif
	1020	! Temperature
	1021	if (guide_T) then
[1188]	1022	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
	1023	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[1170]	1024	print*,'ncidT,varidT',ncidt,varidt
	1025	if (ncidpl.eq.-99) ncidpl=ncidt
	1026	endif
	1027	! Humidite
	1028	if (guide_Q) then
[1188]	1029	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
	1030	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[1170]	1031	print*,'ncidQ,varidQ',ncidQ,varidQ
	1032	if (ncidpl.eq.-99) ncidpl=ncidQ
	1033	endif
	1034	! Pression de surface
	1035	if ((guide_P).OR.(guide_modele)) then
[1188]	1036	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
	1037	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[1170]	1038	print*,'ncidps,varidps',ncidps,varidps
	1039	endif
	1040	! Coordonnee verticale
	1041	if (.not.guide_modele) then
[1188]	1042	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
	1043	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[1170]	1044	print*,'ncidpl,varidpl',ncidpl,varidpl
	1045	endif
	1046	! Coefs ap, bp pour calcul de la pression aux differents niveaux
	1047	if (guide_modele) then
	1048	#ifdef NC_DOUBLE
	1049	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	1050	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	1051	#else
	1052	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	1053	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	1054	#endif
	1055	else
	1056	#ifdef NC_DOUBLE
	1057	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	1058	#else
	1059	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	1060	#endif
	1061	apnc=apnc*100.! conversion en Pascals
	1062	bpnc(:)=0.
	1063	endif
	1064	first=.FALSE.
	1065	endif ! (first)
	1066
	1067	! -----------------------------------------------------------------
	1068	! lecture des champs u, v, T, Q, ps
	1069	! -----------------------------------------------------------------
	1070
	1071	! dimensions pour les champs scalaires et le vent zonal
	1072	start(1)=1
	1073	start(2)=1
	1074	start(3)=1
	1075	start(4)=timestep
	1076
	1077	count(1)=iip1
	1078	count(2)=jjp1
	1079	count(3)=nlevnc
	1080	count(4)=1
	1081
	1082	! Vent zonal
	1083	if (guide_u) then
	1084	#ifdef NC_DOUBLE
	1085	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
	1086	#else
	1087	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
	1088	#endif
	1089	IF (invert_y) THEN
[1304]	1090	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1170]	1091	ENDIF
	1092	endif
	1093
	1094	! Temperature
	1095	if (guide_T) then
	1096	#ifdef NC_DOUBLE
	1097	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
	1098	#else
	1099	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
	1100	#endif
	1101	IF (invert_y) THEN
[1304]	1102	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1170]	1103	ENDIF
	1104	endif
	1105
	1106	! Humidite
	1107	if (guide_Q) then
	1108	#ifdef NC_DOUBLE
	1109	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
	1110	#else
	1111	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
	1112	#endif
	1113	IF (invert_y) THEN
[1304]	1114	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1170]	1115	ENDIF
	1116
	1117	endif
	1118
	1119	! Vent meridien
	1120	if (guide_v) then
	1121	count(2)=jjm
	1122	#ifdef NC_DOUBLE
	1123	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
	1124	#else
	1125	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
	1126	#endif
	1127	IF (invert_y) THEN
[1304]	1128	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1170]	1129	ENDIF
	1130	endif
	1131
	1132	! Pression de surface
	1133	if ((guide_P).OR.(guide_modele)) then
	1134	start(3)=timestep
	1135	start(4)=0
	1136	count(2)=jjp1
	1137	count(3)=1
	1138	count(4)=0
	1139	#ifdef NC_DOUBLE
	1140	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
	1141	#else
	1142	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
	1143	#endif
	1144	IF (invert_y) THEN
	1145	CALL invert_lat(iip1,jjp1,1,psnat2)
	1146	ENDIF
	1147	endif
	1148
	1149	END SUBROUTINE guide_read
	1150
	1151	!=======================================================================
	1152	SUBROUTINE guide_read2D(timestep)
	1153
	1154	IMPLICIT NONE
	1155
	1156	#include "netcdf.inc"
	1157	#include "dimensions.h"
	1158	#include "paramet.h"
	1159
	1160	INTEGER, INTENT(IN) :: timestep
	1161
	1162	LOGICAL, SAVE :: first=.TRUE.
	1163	! Identification fichiers et variables NetCDF:
	1164	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
	1165	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
	1166	INTEGER :: ncidpl,varidpl,varidap,varidbp
	1167	! Variables auxiliaires NetCDF:
	1168	INTEGER, DIMENSION(4) :: start,count
	1169	INTEGER :: status,rcode
	1170	! Variables for 3D extension:
	1171	REAL, DIMENSION (jjp1,llm) :: zu
	1172	REAL, DIMENSION (jjm,llm) :: zv
	1173	INTEGER :: i
	1174
	1175	! -----------------------------------------------------------------
	1176	! Premier appel: initialisation de la lecture des fichiers
	1177	! -----------------------------------------------------------------
	1178	if (first) then
	1179	ncidpl=-99
	1180	print*,'Guide: ouverture des fichiers guidage '
	1181	! Niveaux de pression si non constants
	1182	if (guide_modele) then
	1183	print *,'Lecture du guidage sur niveaux mod�le'
[1188]	1184	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
	1185	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
	1186	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[1170]	1187	print*,'ncidpl,varidap',ncidpl,varidap
	1188	endif
	1189	! Vent zonal
	1190	if (guide_u) then
[1188]	1191	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
	1192	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[1170]	1193	print*,'ncidu,varidu',ncidu,varidu
	1194	if (ncidpl.eq.-99) ncidpl=ncidu
	1195	endif
	1196	! Vent meridien
	1197	if (guide_v) then
[1188]	1198	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
	1199	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[1170]	1200	print*,'ncidv,varidv',ncidv,varidv
	1201	if (ncidpl.eq.-99) ncidpl=ncidv
	1202	endif
	1203	! Temperature
	1204	if (guide_T) then
[1188]	1205	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
	1206	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[1170]	1207	print*,'ncidT,varidT',ncidt,varidt
	1208	if (ncidpl.eq.-99) ncidpl=ncidt
	1209	endif
	1210	! Humidite
	1211	if (guide_Q) then
[1188]	1212	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
	1213	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[1170]	1214	print*,'ncidQ,varidQ',ncidQ,varidQ
	1215	if (ncidpl.eq.-99) ncidpl=ncidQ
	1216	endif
	1217	! Pression de surface
	1218	if ((guide_P).OR.(guide_modele)) then
[1188]	1219	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
	1220	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[1170]	1221	print*,'ncidps,varidps',ncidps,varidps
	1222	endif
	1223	! Coordonnee verticale
	1224	if (.not.guide_modele) then
[1188]	1225	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
	1226	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[1170]	1227	print*,'ncidpl,varidpl',ncidpl,varidpl
	1228	endif
	1229	! Coefs ap, bp pour calcul de la pression aux differents niveaux
	1230	if (guide_modele) then
	1231	#ifdef NC_DOUBLE
	1232	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	1233	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	1234	#else
	1235	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	1236	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	1237	#endif
	1238	else
	1239	#ifdef NC_DOUBLE
	1240	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	1241	#else
	1242	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	1243	#endif
	1244	apnc=apnc*100.! conversion en Pascals
	1245	bpnc(:)=0.
	1246	endif
	1247	first=.FALSE.
	1248	endif ! (first)
	1249
	1250	! -----------------------------------------------------------------
	1251	! lecture des champs u, v, T, Q, ps
	1252	! -----------------------------------------------------------------
	1253
	1254	! dimensions pour les champs scalaires et le vent zonal
	1255	start(1)=1
	1256	start(2)=1
	1257	start(3)=1
	1258	start(4)=timestep
	1259
	1260	count(1)=1
	1261	count(2)=jjp1
	1262	count(3)=nlevnc
	1263	count(4)=1
	1264
	1265	! Vent zonal
	1266	if (guide_u) then
	1267	#ifdef NC_DOUBLE
	1268	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
	1269	#else
	1270	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
	1271	#endif
	1272	DO i=1,iip1
	1273	unat2(i,:,:)=zu(:,:)
	1274	ENDDO
	1275
	1276	IF (invert_y) THEN
[1304]	1277	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1170]	1278	ENDIF
	1279
	1280	endif
	1281
	1282	! Temperature
	1283	if (guide_T) then
	1284	#ifdef NC_DOUBLE
	1285	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
	1286	#else
	1287	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
	1288	#endif
	1289	DO i=1,iip1
	1290	tnat2(i,:,:)=zu(:,:)
	1291	ENDDO
	1292
	1293	IF (invert_y) THEN
[1304]	1294	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1170]	1295	ENDIF
	1296
	1297	endif
	1298
	1299	! Humidite
	1300	if (guide_Q) then
	1301	#ifdef NC_DOUBLE
	1302	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
	1303	#else
	1304	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
	1305	#endif
	1306	DO i=1,iip1
	1307	qnat2(i,:,:)=zu(:,:)
	1308	ENDDO
	1309
	1310	IF (invert_y) THEN
[1304]	1311	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1170]	1312	ENDIF
	1313
	1314	endif
	1315
	1316	! Vent meridien
	1317	if (guide_v) then
	1318	count(2)=jjm
	1319	#ifdef NC_DOUBLE
	1320	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
	1321	#else
	1322	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
	1323	#endif
	1324	DO i=1,iip1
	1325	vnat2(i,:,:)=zv(:,:)
	1326	ENDDO
	1327
	1328	IF (invert_y) THEN
[1304]	1329	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1170]	1330	ENDIF
	1331
	1332	endif
	1333
	1334	! Pression de surface
	1335	if ((guide_P).OR.(guide_modele)) then
	1336	start(3)=timestep
	1337	start(4)=0
	1338	count(2)=jjp1
	1339	count(3)=1
	1340	count(4)=0
	1341	#ifdef NC_DOUBLE
	1342	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
	1343	#else
	1344	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
	1345	#endif
	1346	DO i=1,iip1
	1347	psnat2(i,:)=zu(:,1)
	1348	ENDDO
	1349
	1350	IF (invert_y) THEN
	1351	CALL invert_lat(iip1,jjp1,1,psnat2)
	1352	ENDIF
	1353
	1354	endif
	1355
	1356	END SUBROUTINE guide_read2D
	1357
	1358	!=======================================================================
	1359	SUBROUTINE guide_out(varname,hsize,vsize,field)
	1360
	1361	IMPLICIT NONE
	1362
	1363	INCLUDE "dimensions.h"
	1364	INCLUDE "paramet.h"
	1365	INCLUDE "netcdf.inc"
	1366	INCLUDE "comgeom2.h"
	1367	INCLUDE "comconst.h"
	1368	INCLUDE "comvert.h"
	1369
	1370	! Variables entree
	1371	CHARACTER, INTENT(IN) :: varname
	1372	INTEGER, INTENT (IN) :: hsize,vsize
	1373	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
	1374
	1375	! Variables locales
	1376	INTEGER, SAVE :: timestep=0
	1377	! Identites fichier netcdf
	1378	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
	1379	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
	1380	INTEGER, DIMENSION (3) :: dim3
	1381	INTEGER, DIMENSION (4) :: dim4,count,start
	1382	INTEGER :: ierr, varid
	1383
	1384	print *,'Guide: output timestep',timestep,'var ',varname
	1385	IF (timestep.EQ.0) THEN
	1386	! ----------------------------------------------
	1387	! initialisation fichier de sortie
	1388	! ----------------------------------------------
	1389	! Ouverture du fichier
	1390	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
	1391	! Definition des dimensions
	1392	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
	1393	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
	1394	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
	1395	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
	1396	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
	1397	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
	1398
	1399	! Creation des variables dimensions
	1400	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
	1401	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
	1402	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
	1403	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
	1404	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
	1405	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
	1406	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
	1407
	1408	ierr=NF_ENDDEF(nid)
	1409
	1410	! Enregistrement des variables dimensions
	1411	#ifdef NC_DOUBLE
	1412	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
	1413	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
	1414	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
	1415	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
	1416	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
	1417	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
	1418	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
	1419	#else
	1420	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
	1421	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
	1422	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
	1423	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
	1424	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
	1425	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
	1426	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
	1427	#endif
	1428	! --------------------------------------------------------------------
	1429	! Cr�ation des variables sauvegard�es
	1430	! --------------------------------------------------------------------
	1431	ierr = NF_REDEF(nid)
	1432	! Surface pressure (GCM)
	1433	dim3=(/id_lonv,id_latu,id_tim/)
	1434	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
	1435	! Surface pressure (guidage)
	1436	IF (guide_P) THEN
	1437	dim3=(/id_lonv,id_latu,id_tim/)
	1438	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
	1439	ENDIF
	1440	! Zonal wind
	1441	IF (guide_u) THEN
	1442	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
	1443	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
	1444	ENDIF
	1445	! Merid. wind
	1446	IF (guide_v) THEN
	1447	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
	1448	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
	1449	ENDIF
	1450	! Pot. Temperature
	1451	IF (guide_T) THEN
	1452	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	1453	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
	1454	ENDIF
	1455	! Specific Humidity
	1456	IF (guide_Q) THEN
	1457	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	1458	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
	1459	ENDIF
	1460
	1461	ierr = NF_ENDDEF(nid)
	1462	ierr = NF_CLOSE(nid)
	1463	ENDIF ! timestep=0
	1464
	1465	! --------------------------------------------------------------------
	1466	! Enregistrement du champ
	1467	! --------------------------------------------------------------------
	1468	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
	1469
	1470	SELECT CASE (varname)
	1471	CASE ("S")
	1472	timestep=timestep+1
	1473	ierr = NF_INQ_VARID(nid,"SP",varid)
	1474	start=(/1,1,timestep,0/)
	1475	count=(/iip1,jjp1,1,0/)
	1476	#ifdef NC_DOUBLE
	1477	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1478	#else
	1479	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1480	#endif
	1481	CASE ("P")
	1482	ierr = NF_INQ_VARID(nid,"ps",varid)
	1483	start=(/1,1,timestep,0/)
	1484	count=(/iip1,jjp1,1,0/)
	1485	#ifdef NC_DOUBLE
	1486	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1487	#else
	1488	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1489	#endif
	1490	CASE ("U")
	1491	ierr = NF_INQ_VARID(nid,"ucov",varid)
	1492	start=(/1,1,1,timestep/)
	1493	count=(/iip1,jjp1,llm,1/)
	1494	#ifdef NC_DOUBLE
	1495	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1496	#else
	1497	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1498	#endif
	1499	CASE ("V")
	1500	ierr = NF_INQ_VARID(nid,"vcov",varid)
	1501	start=(/1,1,1,timestep/)
	1502	count=(/iip1,jjm,llm,1/)
	1503	#ifdef NC_DOUBLE
	1504	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1505	#else
	1506	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1507	#endif
	1508	CASE ("T")
	1509	ierr = NF_INQ_VARID(nid,"teta",varid)
	1510	start=(/1,1,1,timestep/)
	1511	count=(/iip1,jjp1,llm,1/)
	1512	#ifdef NC_DOUBLE
	1513	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1514	#else
	1515	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1516	#endif
	1517	CASE ("Q")
	1518	ierr = NF_INQ_VARID(nid,"q",varid)
	1519	start=(/1,1,1,timestep/)
	1520	count=(/iip1,jjp1,llm,1/)
	1521	#ifdef NC_DOUBLE
	1522	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
	1523	#else
	1524	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
	1525	#endif
	1526	END SELECT
	1527
	1528	ierr = NF_CLOSE(nid)
	1529
	1530	END SUBROUTINE guide_out
	1531
	1532
	1533	!===========================================================================
	1534	subroutine correctbid(iim,nl,x)
	1535	integer iim,nl
	1536	real x(iim+1,nl)
	1537	integer i,l
	1538	real zz
	1539
	1540	do l=1,nl
	1541	do i=2,iim-1
	1542	if(abs(x(i,l)).gt.1.e10) then
	1543	zz=0.5*(x(i-1,l)+x(i+1,l))
	1544	print*,'correction ',i,l,x(i,l),zz
	1545	x(i,l)=zz
	1546	endif
	1547	enddo
	1548	enddo
	1549	return
	1550	end subroutine correctbid
	1551
	1552	!===========================================================================
	1553	END MODULE guide_mod

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

Download in other formats:

Original Format