Context Navigation

guide_mod.F90 @ 3020

Last change on this file since 3020 was 2572, checked in by emillour, 3 years ago
Common dynamics: Fixes for the picky gfortran10 compiler which identifies using a scalar instead of a one-element array as an error. MW+EM
File size: 59.3 KB

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

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

Download in other formats:

Original Format