Context Navigation

guide_mod.F90 @ 5411

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

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

Download in other formats:

Original Format