Context Navigation

guide_mod.F90 @ 5183

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

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

Download in other formats:

Original Format