Context Navigation

guide_mod.F90 @ 5155

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

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

Download in other formats:

Original Format