Context Navigation

guide_mod.F90 @ 1678

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

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

Download in other formats:

Original Format