Context Navigation

guide_p_mod.F90 @ 1563

Last change on this file since 1563 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:executable set to *

File size: 64.1 KB

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

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

Download in other formats:

Original Format