Context Navigation

guide_mod.f90 @ 5452

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

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

Download in other formats:

Original Format