Context Navigation

guide_mod.F90 @ 2280

Last change on this file since 2280 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``*
File size: 58.4 KB

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

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

Download in other formats:

Original Format