Context Navigation

guide_p_mod.F90 @ 2544

Last change on this file since 2544 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 71.9 KB

Rev	Line
[1172]	1	!
[1403]	2	! $Id: guide_p_mod.F90 2160 2014-11-28 15:36:29Z fairhead $
[1172]	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
[1750]	13	USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
	14	USE pres2lev_mod
[1172]	15
	16	IMPLICIT NONE
	17
	18	! ---------------------------------------------
	19	! Declarations des cles logiques et parametres
	20	! ---------------------------------------------
	21	INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
[1512]	22	INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
[1172]	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
[1512]	26	LOGICAL, PRIVATE, SAVE :: invert_p,invert_y,ini_anal
	27	LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav,guide_modele
[1172]	28
	29	REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
	30	REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
	31	REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
	32	REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
	33	REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
	34
	35	REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
	36	REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
	37	REAL, PRIVATE, SAVE :: tau_lon,tau_lat
	38
	39	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
	40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q
	41	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
	42
	43	! ---------------------------------------------
	44	! Variables de guidage
	45	! ---------------------------------------------
	46	! Variables des fichiers de guidage
	47	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
	48	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
	49	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
	50	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
[1512]	51	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: pnat1,pnat2
[1172]	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
[2160]	69
[1403]	70	USE control_mod
[2160]	71
[1172]	72	IMPLICIT NONE
	73
	74	INCLUDE "dimensions.h"
	75	INCLUDE "paramet.h"
	76	INCLUDE "netcdf.inc"
	77
[2160]	78	! For grossismx:
	79	include "serre.h"
	80
[1172]	81	INTEGER :: error,ncidpl,rid,rcod
	82	CHARACTER (len = 80) :: abort_message
	83	CHARACTER (len = 20) :: modname = 'guide_init'
	84
	85	! ---------------------------------------------
	86	! Lecture des parametres:
	87	! ---------------------------------------------
	88	! Variables guidees
	89	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
	90	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
	91	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
	92	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
	93	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
	94	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
	95	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
	96
	97	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
	98	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
[2160]	99	if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
	100	call abort_gcm("guide_init", &
	101	"zonal nudging requires grid regular in longitude", 1)
[1172]	102
	103	! Constantes de rappel. Unite : fraction de jour
	104	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
	105	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
	106	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
	107	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
	108	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
	109	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
	110	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
	111	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
	112	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
	113	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
	114	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
	115	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
	116
	117	! Sauvegarde du for�age
	118	CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
	119	CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
	120	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
	121	IF (iguide_sav.GT.0) THEN
[2160]	122	iguide_sav=day_step/iguide_sav
	123	ELSE if (iguide_sav == 0) then
	124	iguide_sav = huge(0)
[1172]	125	ELSE
[2160]	126	iguide_sav=day_step*iguide_sav
[1172]	127	ENDIF
	128
	129	! Guidage regional seulement (sinon constant ou suivant le zoom)
	130	CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
	131	CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
	132	CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
	133	CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
	134	CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
	135	CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
	136	CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
	137
	138	! Parametres pour lecture des fichiers
	139	CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
[1512]	140	CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
	141	IF (iguide_int.EQ.0) THEN
	142	iguide_int=1
	143	ELSEIF (iguide_int.GT.0) THEN
[1172]	144	iguide_int=day_step/iguide_int
	145	ELSE
	146	iguide_int=day_step*iguide_int
	147	ENDIF
[1512]	148	CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
	149	! Pour compatibilite avec ancienne version avec guide_modele
	150	CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
	151	IF (guide_modele) THEN
	152	guide_plevs=1
	153	ENDIF
	154	! Fin raccord
[1172]	155	CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
	156	CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
	157	CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
	158	CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
	159
	160	! ---------------------------------------------
	161	! Determination du nombre de niveaux verticaux
	162	! des fichiers guidage
	163	! ---------------------------------------------
	164	ncidpl=-99
[1512]	165	if (guide_plevs.EQ.1) then
[2160]	166	if (ncidpl.eq.-99) then
[1910]	167	rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
	168	if (rcod.NE.NF_NOERR) THEN
	169	print *,'Guide: probleme -> pas de fichier apbp.nc'
	170	CALL abort_gcm(modname,abort_message,1)
	171	endif
	172	endif
[1512]	173	elseif (guide_plevs.EQ.2) then
[2160]	174	if (ncidpl.EQ.-99) then
[1910]	175	rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
	176	if (rcod.NE.NF_NOERR) THEN
	177	print *,'Guide: probleme -> pas de fichier P.nc'
	178	CALL abort_gcm(modname,abort_message,1)
	179	endif
	180	endif
[1512]	181	elseif (guide_u) then
[1910]	182	if (ncidpl.eq.-99) then
	183	rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
	184	if (rcod.NE.NF_NOERR) THEN
	185	print *,'Guide: probleme -> pas de fichier u.nc'
	186	CALL abort_gcm(modname,abort_message,1)
	187	endif
	188	endif
[1512]	189	elseif (guide_v) then
[1910]	190	if (ncidpl.eq.-99) then
	191	rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
	192	if (rcod.NE.NF_NOERR) THEN
	193	print *,'Guide: probleme -> pas de fichier v.nc'
	194	CALL abort_gcm(modname,abort_message,1)
	195	endif
	196	endif
[1512]	197	elseif (guide_T) then
[1910]	198	if (ncidpl.eq.-99) then
	199	rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
	200	if (rcod.NE.NF_NOERR) THEN
	201	print *,'Guide: probleme -> pas de fichier T.nc'
	202	CALL abort_gcm(modname,abort_message,1)
	203	endif
	204	endif
[1512]	205	elseif (guide_Q) then
[1910]	206	if (ncidpl.eq.-99) then
	207	rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
	208	if (rcod.NE.NF_NOERR) THEN
	209	print *,'Guide: probleme -> pas de fichier hur.nc'
	210	CALL abort_gcm(modname,abort_message,1)
	211	endif
	212	endif
[1172]	213	endif
	214	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
	215	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
	216	IF (error.NE.NF_NOERR) THEN
	217	print *,'Guide: probleme lecture niveaux pression'
	218	CALL abort_gcm(modname,abort_message,1)
	219	ENDIF
	220	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
	221	print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
[1188]	222	rcod = nf90_close(ncidpl)
[1172]	223
	224	! ---------------------------------------------
	225	! Allocation des variables
	226	! ---------------------------------------------
	227	abort_message='pb in allocation guide'
	228
	229	ALLOCATE(apnc(nlevnc), stat = error)
	230	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	231	ALLOCATE(bpnc(nlevnc), stat = error)
	232	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	233	apnc=0.;bpnc=0.
	234
	235	ALLOCATE(alpha_pcor(llm), stat = error)
	236	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	237	ALLOCATE(alpha_u(ip1jmp1), stat = error)
	238	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	239	ALLOCATE(alpha_v(ip1jm), stat = error)
	240	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	241	ALLOCATE(alpha_T(ip1jmp1), stat = error)
	242	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	243	ALLOCATE(alpha_Q(ip1jmp1), stat = error)
	244	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	245	ALLOCATE(alpha_P(ip1jmp1), stat = error)
	246	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	247	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
	248
	249	IF (guide_u) THEN
	250	ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
	251	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	252	ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
	253	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	254	ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
	255	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	256	ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
	257	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	258	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
	259	ENDIF
	260
	261	IF (guide_T) THEN
	262	ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
	263	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	264	ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
	265	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	266	ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
	267	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	268	ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
	269	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	270	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
	271	ENDIF
	272
	273	IF (guide_Q) THEN
	274	ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
	275	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	276	ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
	277	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	278	ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
	279	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	280	ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
	281	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	282	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
	283	ENDIF
	284
	285	IF (guide_v) THEN
	286	ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
	287	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	288	ALLOCATE(vgui1(ip1jm,llm), stat = error)
	289	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	290	ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
	291	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	292	ALLOCATE(vgui2(ip1jm,llm), stat = error)
	293	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	294	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
	295	ENDIF
	296
[1512]	297	IF (guide_plevs.EQ.2) THEN
	298	ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
	299	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	300	ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
	301	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	302	pnat1=0.;pnat2=0.;
	303	ENDIF
	304
	305	IF (guide_P.OR.guide_plevs.EQ.1) THEN
[1172]	306	ALLOCATE(psnat1(iip1,jjp1), stat = error)
	307	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	308	ALLOCATE(psnat2(iip1,jjp1), stat = error)
	309	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	310	psnat1=0.;psnat2=0.;
	311	ENDIF
	312	IF (guide_P) THEN
	313	ALLOCATE(psgui2(ip1jmp1), stat = error)
	314	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	315	ALLOCATE(psgui1(ip1jmp1), stat = error)
	316	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	317	psgui1=0.;psgui2=0.
	318	ENDIF
	319
	320	! ---------------------------------------------
	321	! Lecture du premier etat de guidage.
	322	! ---------------------------------------------
	323	IF (guide_2D) THEN
	324	CALL guide_read2D(1)
	325	ELSE
	326	CALL guide_read(1)
	327	ENDIF
	328	IF (guide_v) vnat1=vnat2
	329	IF (guide_u) unat1=unat2
	330	IF (guide_T) tnat1=tnat2
	331	IF (guide_Q) qnat1=qnat2
[1512]	332	IF (guide_plevs.EQ.2) pnat1=pnat2
	333	IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
[1172]	334
	335	END SUBROUTINE guide_init
	336
	337	!=======================================================================
	338	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
[2056]	339	use exner_hyb_p_m, only: exner_hyb_p
	340	use exner_milieu_p_m, only: exner_milieu_p
[1864]	341	USE parallel_lmdz
[1403]	342	USE control_mod
[1172]	343
	344	IMPLICIT NONE
	345
	346	INCLUDE "dimensions.h"
	347	INCLUDE "paramet.h"
	348	INCLUDE "comconst.h"
	349	INCLUDE "comvert.h"
	350
	351	! Variables entree
	352	INTEGER, INTENT(IN) :: itau !pas de temps
	353	REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
	354	REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
	355	REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
	356
	357	! Variables locales
	358	LOGICAL, SAVE :: first=.TRUE.
	359	LOGICAL :: f_out ! sortie guidage
	360	REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
[1512]	361	! Variables pour fonction Exner (P milieu couche)
[2056]	362	REAL, DIMENSION (iip1,jjp1,llm) :: pk
[1512]	363	REAL, DIMENSION (iip1,jjp1) :: pks
	364	REAL :: unskap
	365	REAL, DIMENSION (ip1jmp1,llmp1) :: p ! besoin si guide_P
[1172]	366	! Compteurs temps:
	367	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
	368	REAL :: ditau, dday_step
	369	REAL :: tau,reste ! position entre 2 etats de guidage
	370	REAL, SAVE :: factt ! pas de temps en fraction de jour
	371
[1512]	372	INTEGER :: i,j,l
[1172]	373
	374	ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1
	375	jjb_u=jj_begin ; jje_u=jj_end ; jjn_u=jje_u-jjb_u+1
	376	ijb_v=ij_begin ; ije_v=ij_end ; ijn_v=ije_v-ijb_v+1
	377	jjb_v=jj_begin ; jje_v=jj_end ; jjn_v=jje_v-jjb_v+1
	378	IF (pole_sud) THEN
	379	ije_v=ij_end-iip1
	380	jje_v=jj_end-1
	381	ijn_v=ije_v-ijb_v+1
	382	jjn_v=jje_v-jjb_v+1
	383	ENDIF
	384
	385	! CALL AllGather_Field(ucov,ip1jmp1,llm)
	386	! CALL AllGather_Field(vcov,ip1jm,llm)
	387	! CALL AllGather_Field(teta,ip1jmp1,llm)
	388	! CALL AllGather_Field(ps,ip1jmp1,1)
	389	! CALL AllGather_Field(q,ip1jmp1,llm)
	390
	391	!-----------------------------------------------------------------------
	392	! Initialisations au premier passage
	393	!-----------------------------------------------------------------------
	394
	395	IF (first) THEN
	396	first=.FALSE.
	397	CALL guide_init
	398	itau_test=1001
	399	step_rea=1
	400	count_no_rea=0
	401	! Calcul des constantes de rappel
	402	factt=dtvr*iperiod/daysec
	403	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
	404	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
	405	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
	406	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
	407	call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
	408	! correction de rappel dans couche limite
	409	if (guide_BL) then
	410	alpha_pcor(:)=1.
	411	else
	412	do l=1,llm
	413	alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
	414	enddo
	415	endif
	416	! ini_anal: etat initial egal au guidage
	417	IF (ini_anal) THEN
	418	CALL guide_interp(ps,teta)
	419	IF (guide_u) ucov(ijb_u:ije_u,:)=ugui2(ijb_u:ije_u,:)
	420	IF (guide_v) vcov(ijb_v:ije_v,:)=ugui2(ijb_v:ije_v,:)
	421	IF (guide_T) teta(ijb_u:ije_u,:)=tgui2(ijb_u:ije_u,:)
	422	IF (guide_Q) q(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)
	423	IF (guide_P) THEN
	424	ps(ijb_u:ije_u)=psgui2(ijb_u:ije_u)
	425	CALL pression_p(ip1jmp1,ap,bp,ps,p)
	426	CALL massdair_p(p,masse)
	427	ENDIF
	428	RETURN
	429	ENDIF
	430	! Verification structure guidage
	431	IF (guide_u) THEN
	432	CALL writefield_p('unat',unat1)
	433	CALL writefield_p('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
	434	ENDIF
	435	IF (guide_T) THEN
	436	CALL writefield_p('tnat',tnat1)
	437	CALL writefield_p('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
	438	ENDIF
	439
	440	ENDIF !first
	441
	442	!-----------------------------------------------------------------------
	443	! Lecture des fichiers de guidage ?
	444	!-----------------------------------------------------------------------
	445	IF (iguide_read.NE.0) THEN
	446	ditau=real(itau)
	447	dday_step=real(day_step)
	448	IF (iguide_read.LT.0) THEN
[1512]	449	tau=ditau/dday_step/REAL(iguide_read)
[1172]	450	ELSE
[1512]	451	tau=REAL(iguide_read)*ditau/dday_step
[1172]	452	ENDIF
	453	reste=tau-AINT(tau)
	454	IF (reste.EQ.0.) THEN
	455	IF (itau_test.EQ.itau) THEN
	456	write(,)'deuxieme passage de advreel a itau=',itau
	457	stop
	458	ELSE
[1418]	459	IF (guide_v) vnat1(:,jjb_v:jje_v,:)=vnat2(:,jjb_v:jje_v,:)
	460	IF (guide_u) unat1(:,jjb_u:jje_u,:)=unat2(:,jjb_u:jje_u,:)
	461	IF (guide_T) tnat1(:,jjb_u:jje_u,:)=tnat2(:,jjb_u:jje_u,:)
	462	IF (guide_Q) qnat1(:,jjb_u:jje_u,:)=qnat2(:,jjb_u:jje_u,:)
[1512]	463	IF (guide_plevs.EQ.2) pnat1(:,jjb_u:jje_u,:)=pnat2(:,jjb_u:jje_u,:)
	464	IF (guide_P.OR.guide_plevs.EQ.1) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u)
[1172]	465	step_rea=step_rea+1
	466	itau_test=itau
	467	print*,'Lecture fichiers guidage, pas ',step_rea, &
	468	'apres ',count_no_rea,' non lectures'
	469	IF (guide_2D) THEN
	470	CALL guide_read2D(step_rea)
	471	ELSE
	472	CALL guide_read(step_rea)
	473	ENDIF
	474	count_no_rea=0
	475	ENDIF
	476	ELSE
	477	count_no_rea=count_no_rea+1
	478
	479	ENDIF
	480	ENDIF !iguide_read=0
	481
	482	!-----------------------------------------------------------------------
	483	! Interpolation et conversion des champs de guidage
	484	!-----------------------------------------------------------------------
	485	IF (MOD(itau,iguide_int).EQ.0) THEN
	486	CALL guide_interp(ps,teta)
	487	ENDIF
	488	! Repartition entre 2 etats de guidage
	489	IF (iguide_read.NE.0) THEN
	490	tau=reste
	491	ELSE
	492	tau=1.
	493	ENDIF
	494
	495	!-----------------------------------------------------------------------
	496	! Ajout des champs de guidage
	497	!-----------------------------------------------------------------------
	498	! Sauvegarde du guidage?
	499	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
[1512]	500	IF (f_out) THEN
	501	! Calcul niveaux pression milieu de couches
[2024]	502	CALL pression_p( ip1jmp1, ap, bp, ps, p )
	503	if (pressure_exner) then
[2056]	504	CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
[2024]	505	else
[2056]	506	CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
[1520]	507	endif
	508	unskap=1./kappa
[2024]	509	DO l = 1, llm
	510	DO j=jjb_u,jje_u
	511	DO i =1, iip1
	512	p(i+(j-1)iip1,l) = preff ( pk(i,j,l)/cpp) ** unskap
	513	ENDDO
	514	ENDDO
	515	ENDDO
	516	CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
[1512]	517	ENDIF
[1172]	518
	519	if (guide_u) then
	520	if (guide_add) then
	521	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)
	522	else
	523	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)-ucov(ijb_u:ije_u,:)
	524	endif
	525
	526	if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
	527	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
[2024]	528	IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
	529	IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)ugui1(:,:)+tauugui2(:,:),factt)
	530	IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
[1172]	531	ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
	532	endif
	533
	534	if (guide_T) then
	535	if (guide_add) then
	536	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)
	537	else
	538	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)-teta(ijb_u:ije_u,:)
	539	endif
	540	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
	541	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
[2024]	542	IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
[1172]	543	teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
	544	endif
	545
	546	if (guide_P) then
	547	if (guide_add) then
	548	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)
	549	else
	550	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)-ps(ijb_u:ije_u)
	551	endif
	552	if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
	553	CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
[2024]	554	IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
[1172]	555	ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
	556	CALL pression_p(ip1jmp1,ap,bp,ps,p)
	557	CALL massdair_p(p,masse)
	558	endif
	559
	560	if (guide_Q) then
	561	if (guide_add) then
	562	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)
	563	else
	564	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)-q(ijb_u:ije_u,:)
	565	endif
	566	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
	567	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
[2024]	568	IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
[1172]	569	q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
	570	endif
	571
	572	if (guide_v) then
	573	if (guide_add) then
	574	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)
	575	else
	576	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)-vcov(ijb_v:ije_v,:)
	577	endif
	578
	579	if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
	580	CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
[2024]	581	IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
	582	IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)vgui1(:,:)+tauvgui2(:,:),factt)
	583	IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
[1172]	584	vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
	585	endif
	586
	587	END SUBROUTINE guide_main
	588
	589	!=======================================================================
	590	SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
	591	! field1=afield1+alphafield2
	592
	593	IMPLICIT NONE
	594	INCLUDE "dimensions.h"
	595	INCLUDE "paramet.h"
	596
	597	! input variables
	598	INTEGER, INTENT(IN) :: hsize
	599	INTEGER, INTENT(IN) :: vsize
	600	REAL, DIMENSION(hsize), INTENT(IN) :: alpha
	601	REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
	602
	603	! Local variables
	604	INTEGER :: l
	605
	606	IF (hsize==ip1jm) THEN
	607	do l=1,vsize
	608	field(ijb_v:ije_v,l)=alpha(ijb_v:ije_v)field(ijb_v:ije_v,l)alpha_pcor(l)
	609	enddo
	610	ELSE
	611	do l=1,vsize
	612	field(ijb_u:ije_u,l)=alpha(ijb_u:ije_u)field(ijb_u:ije_u,l)alpha_pcor(l)
	613	enddo
	614	ENDIF
	615
	616	END SUBROUTINE guide_addfield
	617
	618	!=======================================================================
	619	SUBROUTINE guide_zonave(typ,hsize,vsize,field)
	620
	621	IMPLICIT NONE
	622
	623	INCLUDE "dimensions.h"
	624	INCLUDE "paramet.h"
	625	INCLUDE "comgeom.h"
	626	INCLUDE "comconst.h"
	627
	628	! input/output variables
	629	INTEGER, INTENT(IN) :: typ
	630	INTEGER, INTENT(IN) :: vsize
	631	INTEGER, INTENT(IN) :: hsize
	632	REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
	633
	634	! Local variables
	635	LOGICAL, SAVE :: first=.TRUE.
	636	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
	637	INTEGER :: i,j,l,ij
	638	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
	639	REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
	640
	641	IF (first) THEN
	642	first=.FALSE.
	643	!Compute domain for averaging
	644	lond=rlonu*180./pi
	645	imin(1)=1;imax(1)=iip1;
	646	imin(2)=1;imax(2)=iip1;
	647	IF (guide_reg) THEN
	648	DO i=1,iim
	649	IF (lond(i).LT.lon_min_g) imin(1)=i
	650	IF (lond(i).LE.lon_max_g) imax(1)=i
	651	ENDDO
	652	lond=rlonv*180./pi
	653	DO i=1,iim
	654	IF (lond(i).LT.lon_min_g) imin(2)=i
	655	IF (lond(i).LE.lon_max_g) imax(2)=i
	656	ENDDO
	657	ENDIF
	658	ENDIF
	659
	660	fieldm=0.
	661
	662	IF (hsize==jjm) THEN
	663	DO l=1,vsize
	664	! Compute zonal average
	665	DO j=jjb_v,jje_v
	666	DO i=imin(typ),imax(typ)
	667	ij=(j-1)*iip1+i
	668	fieldm(j,l)=fieldm(j,l)+field(ij,l)
	669	ENDDO
	670	ENDDO
[1512]	671	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
[1172]	672	! Compute forcing
	673	DO j=jjb_v,jje_v
	674	DO i=1,iip1
	675	ij=(j-1)*iip1+i
	676	field(ij,l)=fieldm(j,l)
	677	ENDDO
	678	ENDDO
	679	ENDDO
	680	ELSE
	681	DO l=1,vsize
	682	! Compute zonal average
	683	DO j=jjb_v,jje_v
	684	DO i=imin(typ),imax(typ)
	685	ij=(j-1)*iip1+i
	686	fieldm(j,l)=fieldm(j,l)+field(ij,l)
	687	ENDDO
	688	ENDDO
[1512]	689	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
[1172]	690	! Compute forcing
	691	DO j=jjb_u,jje_u
	692	DO i=1,iip1
	693	ij=(j-1)*iip1+i
	694	field(ij,l)=fieldm(j,l)
	695	ENDDO
	696	ENDDO
	697	ENDDO
	698	ENDIF
	699
	700	END SUBROUTINE guide_zonave
	701
	702	!=======================================================================
	703	SUBROUTINE guide_interp(psi,teta)
[2056]	704	use exner_hyb_p_m, only: exner_hyb_p
	705	use exner_milieu_p_m, only: exner_milieu_p
[1864]	706	USE parallel_lmdz
[1172]	707	USE mod_hallo
	708	USE Bands
	709	IMPLICIT NONE
	710
	711	include "dimensions.h"
	712	include "paramet.h"
	713	include "comvert.h"
	714	include "comgeom2.h"
	715	include "comconst.h"
	716
	717	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
	718	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
	719
	720	LOGICAL, SAVE :: first=.TRUE.
	721	! Variables pour niveaux pression:
	722	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
	723	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
	724	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
	725	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
	726	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
	727	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
	728	REAL, DIMENSION (iip1,jjm,llm) :: pbary
	729	! Variables pour fonction Exner (P milieu couche)
[2056]	730	REAL, DIMENSION (iip1,jjp1,llm) :: pk
[1172]	731	REAL, DIMENSION (iip1,jjp1) :: pks
[1512]	732	REAL :: unskap
[1172]	733	! Pression de vapeur saturante
	734	REAL, DIMENSION (ip1jmp1,llm) :: qsat
	735	!Variables intermediaires interpolation
	736	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
	737	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
	738
	739	INTEGER :: i,j,l,ij
	740	TYPE(Request) :: Req
	741
	742	print *,'Guide: conversion variables guidage'
	743	! -----------------------------------------------------------------
[1512]	744	! Calcul des niveaux de pression champs guidage (pour T et Q)
[1172]	745	! -----------------------------------------------------------------
[1512]	746	IF (guide_plevs.EQ.0) THEN
	747	DO l=1,nlevnc
	748	DO j=jjb_u,jje_u
	749	DO i=1,iip1
	750	plnc2(i,j,l)=apnc(l)
	751	plnc1(i,j,l)=apnc(l)
	752	ENDDO
	753	ENDDO
	754	ENDDO
	755	ENDIF
[1172]	756
	757	if (first) then
	758	first=.FALSE.
	759	print*,'Guide: verification ordre niveaux verticaux'
	760	print*,'LMDZ :'
	761	do l=1,llm
	762	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
	763	+psi(1,jje_u)*(bp(l)+bp(l+1))/2.
	764	enddo
	765	print*,'Fichiers guidage'
[1512]	766	SELECT CASE (guide_plevs)
	767	CASE (0)
	768	do l=1,nlevnc
	769	print*,'PL(',l,')=',plnc2(1,jjb_u,l)
	770	enddo
	771	CASE (1)
	772	DO l=1,nlevnc
	773	print,'PL(',l,')=',apnc(l)+bpnc(l)psnat2(i,jjb_u)
	774	ENDDO
	775	CASE (2)
	776	do l=1,nlevnc
	777	print*,'PL(',l,')=',pnat2(1,jjb_u,l)
	778	enddo
	779	END SELECT
[1172]	780	print *,'inversion de l''ordre: invert_p=',invert_p
	781	if (guide_u) then
	782	do l=1,nlevnc
	783	print*,'U(',l,')=',unat2(1,jjb_u,l)
	784	enddo
	785	endif
	786	if (guide_T) then
	787	do l=1,nlevnc
	788	print*,'T(',l,')=',tnat2(1,jjb_u,l)
	789	enddo
	790	endif
	791	endif
	792
	793	! -----------------------------------------------------------------
	794	! Calcul niveaux pression modele
	795	! -----------------------------------------------------------------
	796
	797	! .... Calcul de pls , pression au milieu des couches ,en Pascals
[1512]	798	IF (guide_plevs.EQ.1) THEN
	799	DO l=1,llm
[2024]	800	DO j=jjb_u,jje_u
	801	DO i =1, iip1
[1512]	802	pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
[2024]	803	ENDDO
	804	ENDDO
[1172]	805	ENDDO
[1512]	806	ELSE
[2024]	807	CALL pression_p( ip1jmp1, ap, bp, psi, p )
	808	if (pressure_exner) then
[2056]	809	CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
[1665]	810	else
[2056]	811	CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
[1520]	812	endif
[2024]	813	unskap=1./kappa
	814	DO l = 1, llm
	815	DO j=jjb_u,jje_u
	816	DO i =1, iip1
	817	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
	818	ENDDO
	819	ENDDO
	820	ENDDO
[1512]	821	ENDIF
[1172]	822
	823	! calcul des pressions pour les grilles u et v
	824	do l=1,llm
	825	do j=jjb_u,jje_u
	826	do i=1,iip1
	827	pext(i,j,l)=pls(i,j,l)*aire(i,j)
	828	enddo
	829	enddo
	830	enddo
	831
[1177]	832	CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
[1172]	833	CALL SendRequest(Req)
	834	CALL WaitRequest(Req)
	835
	836	call massbar_p(pext, pbarx, pbary )
	837	do l=1,llm
	838	do j=jjb_u,jje_u
	839	do i=1,iip1
	840	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
	841	plsnc(i,j,l)=pls(i,j,l)
	842	enddo
	843	enddo
	844	enddo
	845	do l=1,llm
	846	do j=jjb_v,jje_v
	847	do i=1,iip1
	848	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
	849	enddo
	850	enddo
	851	enddo
	852
	853	! -----------------------------------------------------------------
[1512]	854	! Interpolation verticale champs guidage sur niveaux modele
[1172]	855	! Conversion en variables gcm (ucov, vcov...)
	856	! -----------------------------------------------------------------
	857	if (guide_P) then
	858	do j=jjb_u,jje_u
	859	do i=1,iim
	860	ij=(j-1)*iip1+i
	861	psgui1(ij)=psnat1(i,j)
	862	psgui2(ij)=psnat2(i,j)
	863	enddo
	864	psgui1(iip1*j)=psnat1(1,j)
	865	psgui2(iip1*j)=psnat2(1,j)
	866	enddo
	867	endif
	868
[1512]	869	IF (guide_T) THEN
	870	! Calcul des nouvelles valeurs des niveaux de pression du guidage
	871	IF (guide_plevs.EQ.1) THEN
	872	DO l=1,nlevnc
	873	DO j=jjb_u,jje_u
	874	DO i=1,iip1
	875	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
	876	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
	877	ENDDO
	878	ENDDO
	879	ENDDO
	880	ELSE IF (guide_plevs.EQ.2) THEN
	881	DO l=1,nlevnc
	882	DO j=jjb_u,jje_u
	883	DO i=1,iip1
	884	plnc2(i,j,l)=pnat2(i,j,l)
	885	plnc1(i,j,l)=pnat1(i,j,l)
	886	ENDDO
	887	ENDDO
	888	ENDDO
	889	ENDIF
[1172]	890
[1512]	891	! Interpolation verticale
[1177]	892	CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
[1512]	893	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
[1177]	894	CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
[1512]	895	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
[1172]	896
[1512]	897	! Conversion en variables GCM
[1172]	898	do l=1,llm
	899	do j=jjb_u,jje_u
	900	IF (guide_teta) THEN
[1512]	901	do i=1,iim
	902	ij=(j-1)*iip1+i
	903	tgui1(ij,l)=zu1(i,j,l)
	904	tgui2(ij,l)=zu2(i,j,l)
	905	enddo
[1172]	906	ELSE
[1512]	907	do i=1,iim
	908	ij=(j-1)*iip1+i
	909	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
	910	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
	911	enddo
[1172]	912	ENDIF
	913	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
	914	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
	915	enddo
	916	do i=1,iip1
	917	tgui1(i,l)=tgui1(1,l)
	918	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
	919	tgui2(i,l)=tgui2(1,l)
	920	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
	921	enddo
	922	enddo
	923	ENDIF
	924
[1512]	925	IF (guide_Q) THEN
	926	! Calcul des nouvelles valeurs des niveaux de pression du guidage
	927	IF (guide_plevs.EQ.1) THEN
	928	DO l=1,nlevnc
	929	DO j=jjb_u,jje_u
	930	DO i=1,iip1
	931	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
	932	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
	933	ENDDO
	934	ENDDO
	935	ENDDO
	936	ELSE IF (guide_plevs.EQ.2) THEN
	937	DO l=1,nlevnc
	938	DO j=jjb_u,jje_u
	939	DO i=1,iip1
	940	plnc2(i,j,l)=pnat2(i,j,l)
	941	plnc1(i,j,l)=pnat1(i,j,l)
	942	ENDDO
	943	ENDDO
	944	ENDDO
	945	ENDIF
[1172]	946
[1512]	947	! Interpolation verticale
[1177]	948	CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
	949	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
	950	CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
	951	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
[1172]	952
[1512]	953	! Conversion en variables GCM
	954	! On suppose qu'on a la bonne variable dans le fichier de guidage:
	955	! Hum.Rel si guide_hr, Hum.Spec. sinon.
[1172]	956	do l=1,llm
[1512]	957	do j=jjb_u,jje_u
[1172]	958	do i=1,iim
	959	ij=(j-1)*iip1+i
	960	qgui1(ij,l)=zu1(i,j,l)
	961	qgui2(ij,l)=zu2(i,j,l)
	962	enddo
	963	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
	964	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
	965	enddo
	966	do i=1,iip1
	967	qgui1(i,l)=qgui1(1,l)
	968	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
	969	qgui2(i,l)=qgui2(1,l)
	970	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
	971	enddo
	972	enddo
	973	IF (guide_hr) THEN
[1177]	974	CALL q_sat(iip1jjn_ullm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp, &
	975	plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
[1172]	976	qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01 !hum. rel. en %
	977	qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01
	978	ENDIF
	979	ENDIF
	980
[1512]	981	IF (guide_u) THEN
	982	! Calcul des nouvelles valeurs des niveaux de pression du guidage
	983	IF (guide_plevs.EQ.1) THEN
	984	DO l=1,nlevnc
	985	DO j=jjb_u,jje_u
	986	DO i=1,iim
	987	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha1p2(i,j) &
	988	& +psnat2(i+1,j)aire(i+1,j)alpha3p4(i+1,j))/aireu(i,j)
	989	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha1p2(i,j) &
	990	& +psnat1(i+1,j)aire(i+1,j)alpha3p4(i+1,j))/aireu(i,j)
	991	ENDDO
	992	plnc2(iip1,j,l)=plnc2(1,j,l)
	993	plnc1(iip1,j,l)=plnc1(1,j,l)
	994	ENDDO
	995	ENDDO
	996	ELSE IF (guide_plevs.EQ.2) THEN
	997	DO l=1,nlevnc
	998	DO j=jjb_u,jje_u
	999	DO i=1,iim
	1000	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha1p2(i,j) &
	1001	& +pnat2(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
	1002	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha1p2(i,j) &
	1003	& +pnat1(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
	1004	ENDDO
	1005	plnc2(iip1,j,l)=plnc2(1,j,l)
	1006	plnc1(iip1,j,l)=plnc1(1,j,l)
	1007	ENDDO
	1008	ENDDO
	1009	ENDIF
	1010
	1011	! Interpolation verticale
	1012	CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
	1013	plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
	1014	CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
	1015	plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
	1016
	1017	! Conversion en variables GCM
	1018	do l=1,llm
	1019	do j=jjb_u,jje_u
	1020	do i=1,iim
	1021	ij=(j-1)*iip1+i
	1022	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
	1023	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
	1024	enddo
	1025	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
	1026	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
	1027	enddo
	1028	do i=1,iip1
	1029	ugui1(i,l)=0.
	1030	ugui1(ip1jm+i,l)=0.
	1031	ugui2(i,l)=0.
	1032	ugui2(ip1jm+i,l)=0.
	1033	enddo
	1034	enddo
	1035	ENDIF
	1036
	1037	IF (guide_v) THEN
	1038	! Calcul des nouvelles valeurs des niveaux de pression du guidage
	1039	IF (guide_plevs.EQ.1) THEN
[2024]	1040	CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
	1041	CALL SendRequest(Req)
	1042	CALL WaitRequest(Req)
	1043	CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
	1044	CALL SendRequest(Req)
	1045	CALL WaitRequest(Req)
[1512]	1046	DO l=1,nlevnc
	1047	DO j=jjb_v,jje_v
	1048	DO i=1,iip1
	1049	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha2p3(i,j) &
	1050	& +psnat2(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
	1051	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha2p3(i,j) &
	1052	& +psnat1(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
	1053	ENDDO
	1054	ENDDO
	1055	ENDDO
	1056	ELSE IF (guide_plevs.EQ.2) THEN
[2024]	1057	CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
	1058	CALL SendRequest(Req)
	1059	CALL WaitRequest(Req)
	1060	CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
	1061	CALL SendRequest(Req)
	1062	CALL WaitRequest(Req)
[1512]	1063	DO l=1,nlevnc
	1064	DO j=jjb_v,jje_v
	1065	DO i=1,iip1
	1066	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha2p3(i,j) &
	1067	& +pnat2(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
	1068	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha2p3(i,j) &
	1069	& +pnat1(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
	1070	ENDDO
	1071	ENDDO
	1072	ENDDO
	1073	ENDIF
	1074	! Interpolation verticale
	1075	CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm, &
	1076	plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
	1077	CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm, &
	1078	plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
	1079	! Conversion en variables GCM
	1080	do l=1,llm
	1081	do j=jjb_v,jje_v
	1082	do i=1,iim
	1083	ij=(j-1)*iip1+i
	1084	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
	1085	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
	1086	enddo
	1087	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
	1088	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
	1089	enddo
	1090	enddo
	1091	ENDIF
	1092
	1093
[1172]	1094	END SUBROUTINE guide_interp
	1095
	1096	!=======================================================================
	1097	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
	1098
	1099	! Calcul des constantes de rappel alpha (=1/tau)
	1100
	1101	implicit none
	1102
	1103	include "dimensions.h"
	1104	include "paramet.h"
	1105	include "comconst.h"
	1106	include "comgeom2.h"
	1107	include "serre.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
[2160]	1261	if (.not. guide_add) alpha = 1. - exp(- alpha)
	1262
[1172]	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:
[1512]	1278	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
	1279	INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
[1172]	1280	INTEGER :: ncidpl,varidpl,varidap,varidbp
	1281	! Variables auxiliaires NetCDF:
	1282	INTEGER, DIMENSION(4) :: start,count
	1283	INTEGER :: status,rcode
	1284
[1910]	1285	CHARACTER (len = 80) :: abort_message
	1286	CHARACTER (len = 20) :: modname = 'guide_read'
[1172]	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 '
[1512]	1293	! Ap et Bp si Niveaux de pression hybrides
	1294	if (guide_plevs.EQ.1) then
	1295	print *,'Lecture du guidage sur niveaux modele'
[1188]	1296	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[1910]	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
[1188]	1301	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[1910]	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
[1188]	1306	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[1910]	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
[1172]	1311	print*,'ncidpl,varidap',ncidpl,varidap
	1312	endif
[1512]	1313	! Pression si guidage sur niveaux P variables
	1314	if (guide_plevs.EQ.2) then
	1315	rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[1910]	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
[1512]	1320	rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[1910]	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
[1512]	1325	print*,'ncidp,varidp',ncidp,varidp
	1326	if (ncidpl.eq.-99) ncidpl=ncidp
	1327	endif
[1172]	1328	! Vent zonal
	1329	if (guide_u) then
[1188]	1330	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[1910]	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
[1188]	1335	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[1910]	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
[1172]	1340	print*,'ncidu,varidu',ncidu,varidu
	1341	if (ncidpl.eq.-99) ncidpl=ncidu
	1342	endif
	1343	! Vent meridien
	1344	if (guide_v) then
[1188]	1345	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[1910]	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
[1188]	1350	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[1910]	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
[1172]	1355	print*,'ncidv,varidv',ncidv,varidv
	1356	if (ncidpl.eq.-99) ncidpl=ncidv
	1357	endif
	1358	! Temperature
	1359	if (guide_T) then
[1188]	1360	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[1910]	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
[1188]	1365	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[1910]	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
[1172]	1370	print*,'ncidT,varidT',ncidt,varidt
	1371	if (ncidpl.eq.-99) ncidpl=ncidt
	1372	endif
	1373	! Humidite
	1374	if (guide_Q) then
[1188]	1375	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[1910]	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
[1188]	1380	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[1910]	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
[1172]	1385	print*,'ncidQ,varidQ',ncidQ,varidQ
	1386	if (ncidpl.eq.-99) ncidpl=ncidQ
	1387	endif
	1388	! Pression de surface
[1512]	1389	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1188]	1390	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[1910]	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
[1188]	1395	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[1910]	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
[1172]	1400	print*,'ncidps,varidps',ncidps,varidps
	1401	endif
	1402	! Coordonnee verticale
[1512]	1403	if (guide_plevs.EQ.0) then
[1188]	1404	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
	1405	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[1172]	1406	print*,'ncidpl,varidpl',ncidpl,varidpl
	1407	endif
	1408	! Coefs ap, bp pour calcul de la pression aux differents niveaux
[1512]	1409	IF (guide_plevs.EQ.1) THEN
[1172]	1410	#ifdef NC_DOUBLE
	1411	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	1412	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	1413	#else
	1414	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	1415	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	1416	#endif
[1512]	1417	ELSEIF (guide_plevs.EQ.0) THEN
[1172]	1418	#ifdef NC_DOUBLE
	1419	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	1420	#else
	1421	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	1422	#endif
	1423	apnc=apnc*100.! conversion en Pascals
	1424	bpnc(:)=0.
[1512]	1425	ENDIF
[1172]	1426	first=.FALSE.
[1512]	1427	ENDIF ! (first)
[1172]	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
[1512]	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
[1172]	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
[1304]	1464	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1172]	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
[1304]	1477	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1172]	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
[1304]	1489	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1172]	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
[1304]	1503	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1172]	1504	ENDIF
	1505	endif
	1506
	1507	! Pression de surface
[1512]	1508	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1172]	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:
[1512]	1539	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
	1540	INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
[1172]	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
[1910]	1550	CHARACTER (len = 80) :: abort_message
	1551	CHARACTER (len = 20) :: modname = 'guide_read2D'
[1172]	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 '
[1512]	1558	! Ap et Bp si niveaux de pression hybrides
	1559	if (guide_plevs.EQ.1) then
[1172]	1560	print *,'Lecture du guidage sur niveaux mod�le'
[1188]	1561	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[1910]	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
[1188]	1566	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[1910]	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
[1188]	1571	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[1910]	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
[1172]	1576	print*,'ncidpl,varidap',ncidpl,varidap
	1577	endif
[1512]	1578	! Pression
	1579	if (guide_plevs.EQ.2) then
	1580	rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[1910]	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
[1512]	1585	rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[1910]	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
[1512]	1590	print*,'ncidp,varidp',ncidp,varidp
	1591	if (ncidpl.eq.-99) ncidpl=ncidp
	1592	endif
[1172]	1593	! Vent zonal
	1594	if (guide_u) then
[1188]	1595	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[1910]	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
[1188]	1600	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[1910]	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
[1172]	1605	print*,'ncidu,varidu',ncidu,varidu
	1606	if (ncidpl.eq.-99) ncidpl=ncidu
	1607	endif
	1608	! Vent meridien
	1609	if (guide_v) then
[1188]	1610	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[1910]	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
[1188]	1615	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[1910]	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
[1172]	1620	print*,'ncidv,varidv',ncidv,varidv
	1621	if (ncidpl.eq.-99) ncidpl=ncidv
	1622	endif
	1623	! Temperature
	1624	if (guide_T) then
[1188]	1625	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[1910]	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
[1188]	1630	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[1910]	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
[1172]	1635	print*,'ncidT,varidT',ncidt,varidt
	1636	if (ncidpl.eq.-99) ncidpl=ncidt
	1637	endif
	1638	! Humidite
	1639	if (guide_Q) then
[1188]	1640	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[1910]	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
[1188]	1645	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[1910]	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
[1172]	1650	print*,'ncidQ,varidQ',ncidQ,varidQ
	1651	if (ncidpl.eq.-99) ncidpl=ncidQ
	1652	endif
	1653	! Pression de surface
[1512]	1654	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1188]	1655	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[1910]	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
[1188]	1660	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[1910]	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
[1172]	1665	print*,'ncidps,varidps',ncidps,varidps
	1666	endif
	1667	! Coordonnee verticale
[1512]	1668	if (guide_plevs.EQ.0) then
[1188]	1669	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
	1670	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[1172]	1671	print*,'ncidpl,varidpl',ncidpl,varidpl
	1672	endif
	1673	! Coefs ap, bp pour calcul de la pression aux differents niveaux
[1512]	1674	if (guide_plevs.EQ.1) then
[1172]	1675	#ifdef NC_DOUBLE
	1676	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	1677	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	1678	#else
	1679	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	1680	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	1681	#endif
[1512]	1682	elseif (guide_plevs.EQ.0) THEN
[1172]	1683	#ifdef NC_DOUBLE
	1684	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	1685	#else
	1686	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	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
[1512]	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
[1172]	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
[1304]	1736	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1172]	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
[1304]	1752	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1172]	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
[1304]	1768	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1172]	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
[1304]	1785	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1172]	1786	ENDIF
	1787	endif
	1788
	1789	! Pression de surface
[1512]	1790	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1172]	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	!=======================================================================
[1512]	1813	SUBROUTINE guide_out(varname,hsize,vsize,field,factt)
[1864]	1814	USE parallel_lmdz
[1172]	1815	IMPLICIT NONE
	1816
	1817	INCLUDE "dimensions.h"
	1818	INCLUDE "paramet.h"
	1819	INCLUDE "netcdf.inc"
	1820	INCLUDE "comgeom2.h"
	1821	INCLUDE "comconst.h"
	1822	INCLUDE "comvert.h"
	1823
	1824	! Variables entree
[2024]	1825	CHARACTER(), INTENT(IN) :: varname
[1172]	1826	INTEGER, INTENT (IN) :: hsize,vsize
	1827	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
[1512]	1828	REAL, INTENT (IN) :: factt
[1172]	1829
	1830	! Variables locales
	1831	INTEGER, SAVE :: timestep=0
	1832	! Identites fichier netcdf
	1833	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
	1834	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
[2024]	1835	INTEGER :: vid_au,vid_av
	1836	INTEGER :: l
[1172]	1837	INTEGER, DIMENSION (3) :: dim3
	1838	INTEGER, DIMENSION (4) :: dim4,count,start
	1839	INTEGER :: ierr, varid
[2024]	1840	REAL, DIMENSION (iip1,hsize,vsize) :: field2
[1172]	1841
	1842	CALL gather_field(field,iip1*hsize,vsize,0)
	1843
	1844	IF (mpi_rank /= 0) RETURN
	1845
	1846	print *,'Guide: output timestep',timestep,'var ',varname
	1847	IF (timestep.EQ.0) THEN
	1848	! ----------------------------------------------
	1849	! initialisation fichier de sortie
	1850	! ----------------------------------------------
	1851	! Ouverture du fichier
	1852	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
	1853	! Definition des dimensions
	1854	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
[2024]	1855	print*,'id_lonu 1 ',id_lonu
[1172]	1856	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
	1857	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
	1858	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
	1859	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
	1860	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
	1861
	1862	! Creation des variables dimensions
	1863	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
[2024]	1864	print*,'id_lonu 2 ',id_lonu
[1172]	1865	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
	1866	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
	1867	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
	1868	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
	1869	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
	1870	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
[2024]	1871	ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
	1872	ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
[1172]	1873
	1874	ierr=NF_ENDDEF(nid)
	1875
	1876	! Enregistrement des variables dimensions
	1877	#ifdef NC_DOUBLE
[2024]	1878	print,'id_lonu DOUBLE ',id_lonu,rlonu180./pi
[1172]	1879	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
	1880	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
	1881	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
	1882	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
	1883	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
	1884	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
	1885	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
[2024]	1886	ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
	1887	ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
[1172]	1888	#else
[2024]	1889	print,'id_lonu 3 ',id_lonu,rlonu180./pi
[1172]	1890	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
	1891	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
	1892	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
	1893	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
	1894	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
	1895	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
	1896	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
[2024]	1897	ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
	1898	ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
[1172]	1899	#endif
	1900	! --------------------------------------------------------------------
	1901	! Cr�ation des variables sauvegard�es
	1902	! --------------------------------------------------------------------
	1903	ierr = NF_REDEF(nid)
[1512]	1904	! Pressure (GCM)
	1905	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[2024]	1906	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
[1172]	1907	! Surface pressure (guidage)
	1908	IF (guide_P) THEN
	1909	dim3=(/id_lonv,id_latu,id_tim/)
	1910	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
	1911	ENDIF
	1912	! Zonal wind
	1913	IF (guide_u) THEN
[2024]	1914	print*,'id_lonu 4 ',id_lonu,varname
[1172]	1915	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
[2024]	1916	ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
	1917	ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
[1172]	1918	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
	1919	ENDIF
	1920	! Merid. wind
	1921	IF (guide_v) THEN
	1922	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
[2024]	1923	ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
	1924	ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
[1172]	1925	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
	1926	ENDIF
	1927	! Pot. Temperature
	1928	IF (guide_T) THEN
	1929	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	1930	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
	1931	ENDIF
	1932	! Specific Humidity
	1933	IF (guide_Q) THEN
	1934	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	1935	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
	1936	ENDIF
	1937
	1938	ierr = NF_ENDDEF(nid)
	1939	ierr = NF_CLOSE(nid)
	1940	ENDIF ! timestep=0
	1941
	1942	! --------------------------------------------------------------------
	1943	! Enregistrement du champ
	1944	! --------------------------------------------------------------------
[1512]	1945
[1172]	1946	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
	1947
[2024]	1948	IF (varname=="SP") timestep=timestep+1
	1949
	1950	IF (varname=="SP") THEN
	1951	print*,'varname=SP=',varname
	1952	ELSE
	1953	print*,'varname diff SP=',varname
	1954	ENDIF
	1955
	1956
	1957	ierr = NF_INQ_VARID(nid,varname,varid)
[1172]	1958	SELECT CASE (varname)
[2024]	1959	CASE ("SP","ps")
[1512]	1960	start=(/1,1,1,timestep/)
	1961	count=(/iip1,jjp1,llm,1/)
[2024]	1962	CASE ("v","va","vcov")
[1172]	1963	start=(/1,1,1,timestep/)
	1964	count=(/iip1,jjm,llm,1/)
[2024]	1965	CASE DEFAULT
[1172]	1966	start=(/1,1,1,timestep/)
	1967	count=(/iip1,jjp1,llm,1/)
[2024]	1968	END SELECT
	1969
	1970	SELECT CASE (varname)
	1971	CASE("u","ua")
	1972	DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
	1973	field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
	1974	CASE("v","va")
	1975	DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
	1976	CASE DEFAULT
	1977	field2=field
	1978	END SELECT
	1979
[1172]	1980	#ifdef NC_DOUBLE
[2024]	1981	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
[1172]	1982	#else
[2024]	1983	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
[1172]	1984	#endif
	1985
	1986	ierr = NF_CLOSE(nid)
	1987
	1988	END SUBROUTINE guide_out
	1989
	1990
	1991	!===========================================================================
	1992	subroutine correctbid(iim,nl,x)
	1993	integer iim,nl
	1994	real x(iim+1,nl)
	1995	integer i,l
	1996	real zz
	1997
	1998	do l=1,nl
	1999	do i=2,iim-1
	2000	if(abs(x(i,l)).gt.1.e10) then
	2001	zz=0.5*(x(i-1,l)+x(i+1,l))
	2002	print*,'correction ',i,l,x(i,l),zz
	2003	x(i,l)=zz
	2004	endif
	2005	enddo
	2006	enddo
	2007	return
	2008	end subroutine correctbid
	2009
	2010	!===========================================================================
	2011	END MODULE guide_p_mod

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

Download in other formats:

Original Format