Context Navigation

guide_p_mod.F90 @ 2438

Last change on this file since 2438 was 2134, checked in by lguez, 10 years ago

In nudging procedures, replaced explicit Euler integration of nudged
fields by exact integration. This does not change anything if
guide_add is true, but it changes the value of alpha if guide_add is
false. We could have taken into account the variation of the nudging
field during a nudging time step. This would be a small correction. We
choose not to take it into account for the time being. Also, we add a
restriction on zonal nudging: we allow it only for a grid which is
regular in longitude. It does not seem to make sense otherwise and the
exact integration would take more programming for an irregular grid.

In the sequential version, copying the parallel versions, set
iguide_int to 1 when the input value is 0.

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 *

File size: 71.9 KB

Rev	Line
[1172]	1	!
[1403]	2	! $Id$
[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
[1716]	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
[2134]	69
[1403]	70	USE control_mod
[2134]	71
[1172]	72	IMPLICIT NONE
	73
	74	INCLUDE "dimensions.h"
	75	INCLUDE "paramet.h"
	76	INCLUDE "netcdf.inc"
	77
[2134]	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')
[2134]	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
[2124]	122	iguide_sav=day_step/iguide_sav
	123	ELSE if (iguide_sav == 0) then
	124	iguide_sav = huge(0)
[1172]	125	ELSE
[2124]	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
[2134]	166	if (ncidpl.eq.-99) then
[1902]	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
[2134]	174	if (ncidpl.EQ.-99) then
[1902]	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
[1902]	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
[1902]	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
[1902]	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
[1902]	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)
[2021]	339	use exner_hyb_p_m, only: exner_hyb_p
	340	use exner_milieu_p_m, only: exner_milieu_p
[1823]	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)
[2021]	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
[2025]	502	CALL pression_p( ip1jmp1, ap, bp, ps, p )
	503	if (pressure_exner) then
[2021]	504	CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
[2025]	505	else
[2021]	506	CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
[1520]	507	endif
	508	unskap=1./kappa
[2025]	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)
[2025]	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)
[2025]	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)
[2025]	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)
[2025]	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)
[2025]	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)
[2021]	704	use exner_hyb_p_m, only: exner_hyb_p
	705	use exner_milieu_p_m, only: exner_milieu_p
[1823]	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)
[2021]	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
[2025]	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.
[2025]	803	ENDDO
	804	ENDDO
[1172]	805	ENDDO
[1512]	806	ELSE
[2025]	807	CALL pression_p( ip1jmp1, ap, bp, psi, p )
	808	if (pressure_exner) then
[2021]	809	CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
[1625]	810	else
[2021]	811	CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
[1520]	812	endif
[2025]	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
[2025]	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
[2025]	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
[2134]	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
[1902]	1285	CHARACTER (len = 80) :: abort_message
[1906]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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
[1902]	1550	CHARACTER (len = 80) :: abort_message
[1906]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1902]	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)
[1823]	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
[2025]	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
[2025]	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
[2025]	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)
[2025]	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)
[2025]	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)
[2025]	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
[2025]	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)
[2025]	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
[2025]	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)
[2025]	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/)
[2025]	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
[2025]	1914	print*,'id_lonu 4 ',id_lonu,varname
[1172]	1915	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
[2025]	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/)
[2025]	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
[2025]	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)
[2025]	1959	CASE ("SP","ps")
[1512]	1960	start=(/1,1,1,timestep/)
	1961	count=(/iip1,jjp1,llm,1/)
[2025]	1962	CASE ("v","va","vcov")
[1172]	1963	start=(/1,1,1,timestep/)
	1964	count=(/iip1,jjm,llm,1/)
[2025]	1965	CASE DEFAULT
[1172]	1966	start=(/1,1,1,timestep/)
	1967	count=(/iip1,jjp1,llm,1/)
[2025]	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
[2025]	1981	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
[1172]	1982	#else
[2025]	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