Context Navigation

guide_p_mod.F90 @ 2902

Last change on this file since 2902 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: 72.5 KB

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

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

Download in other formats:

Original Format