Context Navigation

guide_mod.F90 @ 5496

Last change on this file since 5496 was 5159, checked in by abarral, 6 months ago
Put dimensions.h and paramet.h into 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: 64.3 KB

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

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

Download in other formats:

Original Format