Context Navigation

guide_mod.F90 @ 5095

Last change on this file since 5095 was 5093, checked in by abarral, 4 months ago

Use latest FCM source (2021.05.0) [Note: we still use the legacy FCM1 build system]
Correct UTF8 encoding of french chars

Compil OK (tested: oldrad/rrtm/ecrad, para/seq/1D)
Convergence (ref r5063) bench 33x OK oldrad orch2.0 (tested: para/seq)

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

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

Download in other formats:

Original Format