Context Navigation

guide_mod.F90 @ 5104

Last change on this file since 5104 was 5084, checked in by Laurent Fairhead, 12 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