Context Navigation

guide_mod.F90 @ 4248

Last change on this file since 4248 was 4246, checked in by evignon, 22 months ago
controle dans les .def du presnivs limite en dessous duquel on ne guide plus si guide_BL=false Etienne, pour valentin
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
File size: 66.8 KB

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

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

Download in other formats:

Original Format