Context Navigation

guide_loc_mod.F90 @ 4095

Last change on this file since 4095 was 2740, checked in by lguez, 8 years ago
Added output of nudging coefficients for temperature and humidity, along those for the wind.
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: 82.4 KB

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

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

Download in other formats:

Original Format