Context Navigation

guide_p_mod.F90 @ 1243

Last change on this file since 1243 was 1019, checked in by emillour, 11 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

General stuff:
makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

bld.cfg: adaptations to enable compiling using multiple threads
build_gcm: adaptations to enable compiling using multiple threads
makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change

bibio:
wxios.F90 : Added for possible future use of XIOS library

filtrez:
mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
filtreg_mod : limit use of FFT to parallel mode
mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz

dyn3d:
abort_gcm.F : add things for xios
advtrac.F90 : minor change in CFL outputs
ce0l.F90 : indicesol.h is now module indice_sol_mod
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
inigeom.F : move two computations outside loop

dyn3dpar:
parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
leapfrog_p.F : add xios + correction for times in Newtonian case
ce0l.F90 : indicesol.h is now module indice_sol_mod
inigeom.F : move two computations outside loop

File size: 64.1 KB

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

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

Download in other formats:

Original Format