Context Navigation

guide_p_mod.F90 @ 5002

Last change on this file since 5002 was 1418, checked in by jghattas, 14 years ago

Following bug corrections are already done at branch LMDZ4_AR5 rev
1417,1416 :

Bug corrections for nudged run (pres2lev.F90, guide_p_mod.F90) :

now the results are the same for sequentiel and parallel mode(if adjust=n and compiled with vsafe at mercure).
the results are the same as the sequential mode in previous revision.
test done only with guide_u=y,guide_v=y
copied optimized version of pres2lev.F90 from dyn3dpar to dyn3d

Added condition read_climoz for the variable O3daySTD(calcul_STDlev.h,
undefSTD.F)

Corrected bug in calculation of the diagnostic variable ec550aer
"Extinction at 550nm" (aeropt_5wv.F90) (Maria Raffaella Vuolo, LSCE)

Added stop if use_filtre_fft=y with dyn3d (conf_gcm.F) : this option is
not implemented in dyn3d.

Property svn:executable set to *

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

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

Download in other formats:

Original Format