Context Navigation

guide_loc_mod.F90 @ 5166

Last change on this file since 5166 was 5084, checked in by Laurent Fairhead, 4 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 86.5 KB

Rev	Line
[5084]	1	!
	2	! $Id$
	3	!
[1632]	4	MODULE guide_loc_mod
	5
	6	!=======================================================================
	7	! Auteur: F.Hourdin
	8	! F. Codron 01/09
	9	!=======================================================================
	10
[3995]	11	USE getparam, only: ini_getparam, fin_getparam, getpar
[1632]	12	USE Write_Field_loc
[5084]	13	use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
	14	nf90_inq_dimid, nf90_inquire_dimension
[1823]	15	USE parallel_lmdz
[3995]	16	USE pres2lev_mod, only: pres2lev
[1632]	17
	18	IMPLICIT NONE
	19
	20	! ---------------------------------------------
	21	! Declarations des cles logiques et parametres
	22	! ---------------------------------------------
	23	INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
	24	INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
	25	LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P
	26	LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta
	27	LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
	28	LOGICAL, PRIVATE, SAVE :: invert_p,invert_y,ini_anal
	29	LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav,guide_modele
[3795]	30	!FC
	31	LOGICAL, PRIVATE, SAVE :: convert_Pa
[1632]	32
	33	REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
	34	REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
	35	REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
	36	REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
	37	REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
	38
	39	REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
	40	REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
	41	REAL, PRIVATE, SAVE :: tau_lon,tau_lat
	42
[4246]	43	REAL, PRIVATE, SAVE :: plim_guide_BL
	44
[1632]	45	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
	46	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q
	47	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
	48
	49	! ---------------------------------------------
	50	! Variables de guidage
	51	! ---------------------------------------------
	52	! Variables des fichiers de guidage
	53	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
	54	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
	55	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
	56	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
	57	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: pnat1,pnat2
	58	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2
	59	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc
	60	! Variables aux dimensions du modele
	61	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2
	62	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2
	63	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2
	64	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2
	65	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2
	66
[3995]	67	INTEGER,SAVE,PRIVATE :: ijbu,ijbv,ijeu,ijev !,ijnu,ijnv
[1632]	68	INTEGER,SAVE,PRIVATE :: jjbu,jjbv,jjeu,jjev,jjnu,jjnv
	69
	70
	71	CONTAINS
	72	!=======================================================================
	73
	74	SUBROUTINE guide_init
[2134]	75
[2598]	76	USE control_mod, ONLY: day_step
	77	USE serre_mod, ONLY: grossismx
[2134]	78
[1632]	79	IMPLICIT NONE
	80
	81	INCLUDE "dimensions.h"
	82	INCLUDE "paramet.h"
[5084]	83	INCLUDE "netcdf.inc"
[1632]	84
	85	INTEGER :: error,ncidpl,rid,rcod
	86	CHARACTER (len = 80) :: abort_message
	87	CHARACTER (len = 20) :: modname = 'guide_init'
[3984]	88	CHARACTER (len = 20) :: namedim
[1632]	89
	90	! ---------------------------------------------
	91	! Lecture des parametres:
	92	! ---------------------------------------------
[2263]	93	call ini_getparam("nudging_parameters_out.txt")
[1632]	94	! Variables guidees
	95	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
	96	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
	97	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
	98	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
	99	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
	100	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
	101	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
	102
	103	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
	104	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
[2134]	105	if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
	106	call abort_gcm("guide_init", &
	107	"zonal nudging requires grid regular in longitude", 1)
[1632]	108
	109	! Constantes de rappel. Unite : fraction de jour
	110	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
	111	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
	112	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
	113	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
	114	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
	115	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
	116	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
	117	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
	118	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
	119	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
	120	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
	121	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
[4246]	122	CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
	123
[1632]	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')
	127	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
[5084]	128	IF (iguide_sav.GT.0) THEN
[2124]	129	iguide_sav=day_step/iguide_sav
	130	ELSE if (iguide_sav == 0) then
	131	iguide_sav = huge(0)
[1632]	132	ELSE
[2124]	133	iguide_sav=day_step*iguide_sav
[1632]	134	ENDIF
	135
	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 ')
	144
	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')
[5084]	148	IF (iguide_int.EQ.0) THEN
[1632]	149	iguide_int=1
[5084]	150	ELSEIF (iguide_int.GT.0) THEN
[1632]	151	iguide_int=day_step/iguide_int
	152	ELSE
	153	iguide_int=day_step*iguide_int
	154	ENDIF
	155	CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
	156	! Pour compatibilite avec ancienne version avec guide_modele
	157	CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
	158	IF (guide_modele) THEN
	159	guide_plevs=1
	160	ENDIF
[3795]	161	!FC
	162	CALL getpar('convert_Pa',.true.,convert_Pa,'Convert Pressure levels in Pa')
[1632]	163	! Fin raccord
	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')
	168
[2263]	169	call fin_getparam
	170
[1632]	171	! ---------------------------------------------
	172	! Determination du nombre de niveaux verticaux
	173	! des fichiers guidage
	174	! ---------------------------------------------
	175	ncidpl=-99
[5084]	176	if (guide_plevs.EQ.1) then
	177	if (ncidpl.eq.-99) then
[1902]	178	rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
[5084]	179	if (rcod.NE.NF_NOERR) THEN
[3995]	180	abort_message=' Nudging error -> no file apbp.nc'
[1902]	181	CALL abort_gcm(modname,abort_message,1)
	182	endif
	183	endif
[5084]	184	elseif (guide_plevs.EQ.2) then
	185	if (ncidpl.EQ.-99) then
[1902]	186	rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
[5084]	187	if (rcod.NE.NF_NOERR) THEN
[3995]	188	abort_message=' Nudging error -> no file P.nc'
[1902]	189	CALL abort_gcm(modname,abort_message,1)
	190	endif
	191	endif
[3984]	192
[1632]	193	elseif (guide_u) then
[5084]	194	if (ncidpl.eq.-99) then
[1902]	195	rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
[5084]	196	if (rcod.NE.NF_NOERR) THEN
[3995]	197	abort_message=' Nudging error -> no file u.nc'
[1902]	198	CALL abort_gcm(modname,abort_message,1)
	199	endif
[3984]	200
[1902]	201	endif
[3984]	202
	203
[1632]	204	elseif (guide_v) then
[5084]	205	if (ncidpl.eq.-99) then
[1902]	206	rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
[5084]	207	if (rcod.NE.NF_NOERR) THEN
[3995]	208	abort_message=' Nudging error -> no file v.nc'
[1902]	209	CALL abort_gcm(modname,abort_message,1)
	210	endif
	211	endif
[3984]	212
	213
[1632]	214	elseif (guide_T) then
[5084]	215	if (ncidpl.eq.-99) then
[1902]	216	rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
[5084]	217	if (rcod.NE.NF_NOERR) THEN
[3995]	218	abort_message=' Nudging error -> no file T.nc'
[1902]	219	CALL abort_gcm(modname,abort_message,1)
	220	endif
	221	endif
[3984]	222
	223
	224
[1632]	225	elseif (guide_Q) then
[5084]	226	if (ncidpl.eq.-99) then
[1902]	227	rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
[5084]	228	if (rcod.NE.NF_NOERR) THEN
[3995]	229	abort_message=' Nudging error -> no file hur.nc'
[1902]	230	CALL abort_gcm(modname,abort_message,1)
	231	endif
	232	endif
[3984]	233
	234
[1632]	235	endif
	236	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
[5084]	237	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
	238	IF (error.NE.NF_NOERR) THEN
[3995]	239	abort_message='Nudging: error reading pressure levels'
[1632]	240	CALL abort_gcm(modname,abort_message,1)
	241	ENDIF
	242	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
[3995]	243	write(,)trim(modname)//' : number of vertical levels nlevnc', nlevnc
[1632]	244	rcod = nf90_close(ncidpl)
	245
	246	! ---------------------------------------------
	247	! Allocation des variables
	248	! ---------------------------------------------
[3995]	249	abort_message='nudging allocation error'
[1632]	250
	251	ALLOCATE(apnc(nlevnc), stat = error)
	252	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	253	ALLOCATE(bpnc(nlevnc), stat = error)
	254	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	255	apnc=0.;bpnc=0.
	256
	257	ALLOCATE(alpha_pcor(llm), stat = error)
	258	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	259	ALLOCATE(alpha_u(ijb_u:ije_u), stat = error)
	260	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	261	ALLOCATE(alpha_v(ijb_v:ije_v), stat = error)
	262	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	263	ALLOCATE(alpha_T(ijb_u:ije_u), stat = error)
	264	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	265	ALLOCATE(alpha_Q(ijb_u:ije_u), stat = error)
	266	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	267	ALLOCATE(alpha_P(ijb_u:ije_u), stat = error)
	268	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	269	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
	270
	271	IF (guide_u) THEN
	272	ALLOCATE(unat1(iip1,jjb_u:jje_u,nlevnc), stat = error)
	273	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	274	ALLOCATE(ugui1(ijb_u:ije_u,llm), stat = error)
	275	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	276	ALLOCATE(unat2(iip1,jjb_u:jje_u,nlevnc), stat = error)
	277	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	278	ALLOCATE(ugui2(ijb_u:ije_u,llm), stat = error)
	279	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	280	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
	281	ENDIF
	282
	283	IF (guide_T) THEN
	284	ALLOCATE(tnat1(iip1,jjb_u:jje_u,nlevnc), stat = error)
	285	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	286	ALLOCATE(tgui1(ijb_u:ije_u,llm), stat = error)
	287	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	288	ALLOCATE(tnat2(iip1,jjb_u:jje_u,nlevnc), stat = error)
	289	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	290	ALLOCATE(tgui2(ijb_u:ije_u,llm), stat = error)
	291	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	292	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
	293	ENDIF
	294
	295	IF (guide_Q) THEN
	296	ALLOCATE(qnat1(iip1,jjb_u:jje_u,nlevnc), stat = error)
	297	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	298	ALLOCATE(qgui1(ijb_u:ije_u,llm), stat = error)
	299	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	300	ALLOCATE(qnat2(iip1,jjb_u:jje_u,nlevnc), stat = error)
	301	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	302	ALLOCATE(qgui2(ijb_u:ije_u,llm), stat = error)
	303	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	304	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
	305	ENDIF
	306
	307	IF (guide_v) THEN
	308	ALLOCATE(vnat1(iip1,jjb_v:jje_v,nlevnc), stat = error)
	309	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	310	ALLOCATE(vgui1(ijb_v:ije_v,llm), stat = error)
	311	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	312	ALLOCATE(vnat2(iip1,jjb_v:jje_v,nlevnc), stat = error)
	313	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	314	ALLOCATE(vgui2(ijb_v:ije_v,llm), stat = error)
	315	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	316	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
	317	ENDIF
	318
[5084]	319	IF (guide_plevs.EQ.2) THEN
[1632]	320	ALLOCATE(pnat1(iip1,jjb_u:jje_u,nlevnc), stat = error)
	321	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	322	ALLOCATE(pnat2(iip1,jjb_u:jje_u,nlevnc), stat = error)
	323	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	324	pnat1=0.;pnat2=0.;
	325	ENDIF
	326
[5084]	327	IF (guide_P.OR.guide_plevs.EQ.1) THEN
[1632]	328	ALLOCATE(psnat1(iip1,jjb_u:jje_u), stat = error)
	329	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	330	ALLOCATE(psnat2(iip1,jjb_u:jje_u), stat = error)
	331	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	332	psnat1=0.;psnat2=0.;
	333	ENDIF
	334	IF (guide_P) THEN
	335	ALLOCATE(psgui2(ijb_u:ije_u), stat = error)
	336	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	337	ALLOCATE(psgui1(ijb_u:ije_u), stat = error)
	338	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
	339	psgui1=0.;psgui2=0.
	340	ENDIF
	341
	342	! ---------------------------------------------
	343	! Lecture du premier etat de guidage.
	344	! ---------------------------------------------
	345	IF (guide_2D) THEN
	346	CALL guide_read2D(1)
	347	ELSE
	348	CALL guide_read(1)
	349	ENDIF
	350	IF (guide_v) vnat1=vnat2
	351	IF (guide_u) unat1=unat2
	352	IF (guide_T) tnat1=tnat2
	353	IF (guide_Q) qnat1=qnat2
[5084]	354	IF (guide_plevs.EQ.2) pnat1=pnat2
	355	IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
[1632]	356
	357	END SUBROUTINE guide_init
	358
	359	!=======================================================================
	360	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
[2021]	361	use exner_hyb_loc_m, only: exner_hyb_loc
	362	use exner_milieu_loc_m, only: exner_milieu_loc
[1823]	363	USE parallel_lmdz
[1632]	364	USE control_mod
[1806]	365	USE write_field_loc
[2597]	366	USE comconst_mod, ONLY: cpp, daysec, dtvr, kappa
[2600]	367	USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
[1632]	368
	369	IMPLICIT NONE
	370
	371	INCLUDE "dimensions.h"
	372	INCLUDE "paramet.h"
	373
	374	! Variables entree
	375	INTEGER, INTENT(IN) :: itau !pas de temps
	376	REAL, DIMENSION (ijb_u:ije_u,llm), INTENT(INOUT) :: ucov,teta,q,masse
	377	REAL, DIMENSION (ijb_v:ije_v,llm), INTENT(INOUT) :: vcov
	378	REAL, DIMENSION (ijb_u:ije_u), INTENT(INOUT) :: ps
	379
	380	! Variables locales
	381	LOGICAL, SAVE :: first=.TRUE.
	382	!$OMP THREADPRIVATE(first)
	383	LOGICAL :: f_out ! sortie guidage
[1806]	384	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addu ! var aux: champ de guidage
	385	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addv ! var aux: champ de guidage
[1632]	386	! Variables pour fonction Exner (P milieu couche)
[2021]	387	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pk
[1806]	388	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: pks
[1632]	389	REAL :: unskap
[1806]	390	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: p ! besoin si guide_P
[1632]	391	! Compteurs temps:
	392	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
	393	!$OMP THREADPRIVATE(step_rea,count_no_rea,itau_test)
	394	REAL :: ditau, dday_step
	395	REAL :: tau,reste ! position entre 2 etats de guidage
	396	REAL, SAVE :: factt ! pas de temps en fraction de jour
	397	!$OMP THREADPRIVATE(factt)
	398
	399	INTEGER :: i,j,l
[3995]	400	CHARACTER(LEN=20) :: modname="guide_main"
[1806]	401
[1632]	402	!$OMP MASTER
[3995]	403	ijbu=ij_begin ; ijeu=ij_end
[1632]	404	jjbu=jj_begin ; jjeu=jj_end ; jjnu=jjeu-jjbu+1
[3995]	405	ijbv=ij_begin ; ijev=ij_end
[1632]	406	jjbv=jj_begin ; jjev=jj_end ; jjnv=jjev-jjbv+1
	407	IF (pole_sud) THEN
[3995]	408	ijeu=ij_end-iip1
[1632]	409	ijev=ij_end-iip1
	410	jjev=jj_end-1
	411	jjnv=jjev-jjbv+1
	412	ENDIF
[3995]	413	IF (pole_nord) THEN
	414	ijbu=ij_begin+iip1
	415	ijbv=ij_begin
	416	ENDIF
[1632]	417	!$OMP END MASTER
	418	!$OMP BARRIER
	419
[2036]	420	! PRINT *,'---> on rentre dans guide_main'
[1632]	421	! CALL AllGather_Field(ucov,ip1jmp1,llm)
	422	! CALL AllGather_Field(vcov,ip1jm,llm)
	423	! CALL AllGather_Field(teta,ip1jmp1,llm)
	424	! CALL AllGather_Field(ps,ip1jmp1,1)
	425	! CALL AllGather_Field(q,ip1jmp1,llm)
	426
	427	!-----------------------------------------------------------------------
	428	! Initialisations au premier passage
	429	!-----------------------------------------------------------------------
	430
	431	IF (first) THEN
	432	first=.FALSE.
	433	!$OMP MASTER
[1806]	434	ALLOCATE(f_addu(ijb_u:ije_u,llm) )
	435	ALLOCATE(f_addv(ijb_v:ije_v,llm) )
	436	ALLOCATE(pk(iip1,jjb_u:jje_u,llm) )
	437	ALLOCATE(pks(iip1,jjb_u:jje_u) )
	438	ALLOCATE(p(ijb_u:ije_u,llmp1) )
[1632]	439	CALL guide_init
	440	!$OMP END MASTER
	441	!$OMP BARRIER
	442	itau_test=1001
	443	step_rea=1
	444	count_no_rea=0
	445	! Calcul des constantes de rappel
	446	factt=dtvr*iperiod/daysec
	447	!$OMP MASTER
[1806]	448	call tau2alpha(3, iip1, jjb_v, jje_v, factt, tau_min_v, tau_max_v, alpha_v)
	449	call tau2alpha(2, iip1, jjb_u, jje_u, factt, tau_min_u, tau_max_u, alpha_u)
	450	call tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_T, tau_max_T, alpha_T)
	451	call tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_P, tau_max_P, alpha_P)
	452	call tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_Q, tau_max_Q, alpha_Q)
[1632]	453	! correction de rappel dans couche limite
	454	if (guide_BL) then
	455	alpha_pcor(:)=1.
	456	else
	457	do l=1,llm
[4246]	458	alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
[1632]	459	enddo
	460	endif
	461	!$OMP END MASTER
[1806]	462	!$OMP BARRIER
[1632]	463	! ini_anal: etat initial egal au guidage
	464	IF (ini_anal) THEN
	465	CALL guide_interp(ps,teta)
[2036]	466	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	467	DO l=1,llm
	468	IF (guide_u) ucov(ijbu:ijeu,l)=ugui2(ijbu:ijeu,l)
	469	IF (guide_v) vcov(ijbv:ijev,l)=ugui2(ijbv:ijev,l)
	470	IF (guide_T) teta(ijbu:ijeu,l)=tgui2(ijbu:ijeu,l)
	471	IF (guide_Q) q(ijbu:ijeu,l)=qgui2(ijbu:ijeu,l)
	472	ENDDO
	473
[1632]	474	IF (guide_P) THEN
[1806]	475	!$OMP MASTER
[1632]	476	ps(ijbu:ijeu)=psgui2(ijbu:ijeu)
[1806]	477	!$OMP END MASTER
	478	!$OMP BARRIER
[1632]	479	CALL pression_loc(ijnb_u,ap,bp,ps,p)
	480	CALL massdair_loc(p,masse)
[1806]	481	!$OMP BARRIER
[1632]	482	ENDIF
	483	RETURN
	484	ENDIF
	485
	486	ENDIF !first
	487
	488	!-----------------------------------------------------------------------
	489	! Lecture des fichiers de guidage ?
	490	!-----------------------------------------------------------------------
[5084]	491	IF (iguide_read.NE.0) THEN
[1632]	492	ditau=real(itau)
	493	dday_step=real(day_step)
[5084]	494	IF (iguide_read.LT.0) THEN
[1632]	495	tau=ditau/dday_step/REAL(iguide_read)
	496	ELSE
	497	tau=REAL(iguide_read)*ditau/dday_step
	498	ENDIF
	499	reste=tau-AINT(tau)
[5084]	500	IF (reste.EQ.0.) THEN
	501	IF (itau_test.EQ.itau) THEN
[3995]	502	write(,)trim(modname)//' second pass in advreel at itau=',&
	503	itau
[4469]	504	CALL abort_gcm("guide_loc_lod","stopped",1)
[1632]	505	ELSE
[1806]	506	!$OMP MASTER
[1632]	507	IF (guide_v) vnat1(:,jjbv:jjev,:)=vnat2(:,jjbv:jjev,:)
	508	IF (guide_u) unat1(:,jjbu:jjeu,:)=unat2(:,jjbu:jjeu,:)
	509	IF (guide_T) tnat1(:,jjbu:jjeu,:)=tnat2(:,jjbu:jjeu,:)
	510	IF (guide_Q) qnat1(:,jjbu:jjeu,:)=qnat2(:,jjbu:jjeu,:)
[5084]	511	IF (guide_plevs.EQ.2) pnat1(:,jjbu:jjeu,:)=pnat2(:,jjbu:jjeu,:)
	512	IF (guide_P.OR.guide_plevs.EQ.1) psnat1(:,jjbu:jjeu)=psnat2(:,jjbu:jjeu)
[1806]	513	!$OMP END MASTER
	514	!$OMP BARRIER
[1632]	515	step_rea=step_rea+1
	516	itau_test=itau
[3995]	517	if (is_master) then
	518	write(,)trim(modname)//' Reading nudging files, step ',&
	519	step_rea,'after ',count_no_rea,' skips'
	520	endif
[1632]	521	IF (guide_2D) THEN
[1806]	522	!$OMP MASTER
[1632]	523	CALL guide_read2D(step_rea)
[1806]	524	!$OMP END MASTER
	525	!$OMP BARRIER
[1632]	526	ELSE
[1806]	527	!$OMP MASTER
[1632]	528	CALL guide_read(step_rea)
[1806]	529	!$OMP END MASTER
	530	!$OMP BARRIER
[1632]	531	ENDIF
	532	count_no_rea=0
	533	ENDIF
	534	ELSE
	535	count_no_rea=count_no_rea+1
	536
	537	ENDIF
	538	ENDIF !iguide_read=0
	539
	540	!-----------------------------------------------------------------------
	541	! Interpolation et conversion des champs de guidage
	542	!-----------------------------------------------------------------------
[5084]	543	IF (MOD(itau,iguide_int).EQ.0) THEN
[1632]	544	CALL guide_interp(ps,teta)
	545	ENDIF
	546	! Repartition entre 2 etats de guidage
[5084]	547	IF (iguide_read.NE.0) THEN
[1632]	548	tau=reste
	549	ELSE
	550	tau=1.
	551	ENDIF
	552
[1806]	553	! CALL WriteField_u('ucov_guide',ucov)
	554	! CALL WriteField_v('vcov_guide',vcov)
	555	! CALL WriteField_u('teta_guide',teta)
	556	! CALL WriteField_u('masse_guide',masse)
	557
	558
[3995]	559	!-----------------------------------------------------------------------
[1632]	560	! Ajout des champs de guidage
	561	!-----------------------------------------------------------------------
	562	! Sauvegarde du guidage?
[5084]	563	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
[1632]	564	IF (f_out) THEN
[2036]	565
	566	!$OMP BARRIER
	567	CALL pression_loc(ijnb_u,ap,bp,ps,p)
	568
	569	!$OMP BARRIER
	570	if (pressure_exner) then
	571	CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk)
	572	else
	573	CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk )
	574	endif
	575
	576	!$OMP BARRIER
	577
[1632]	578	unskap=1./kappa
[2036]	579	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2028]	580	DO l = 1, llm
	581	DO j=jjbu,jjeu
	582	DO i =1, iip1
	583	p(i+(j-1)iip1,l) = preff ( pk(i,j,l)/cpp) ** unskap
	584	ENDDO
	585	ENDDO
	586	ENDDO
[2036]	587
	588	CALL guide_out("SP",jjp1,llm,p(ijb_u:ije_u,1:llm),1.)
[1632]	589	ENDIF
	590
	591	if (guide_u) then
	592	if (guide_add) then
[2036]	593	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	594	DO l=1,llm
	595	f_addu(ijbu:ijeu,l)=(1.-tau)ugui1(ijbu:ijeu,l)+tauugui2(ijbu:ijeu,l)
	596	ENDDO
[1632]	597	else
[2036]	598	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	599	DO l=1,llm
	600	f_addu(ijbu:ijeu,l)=(1.-tau)ugui1(ijbu:ijeu,l)+tauugui2(ijbu:ijeu,l)-ucov(ijbu:ijeu,l)
	601	ENDDO
[1632]	602	endif
[1806]	603
	604	! CALL WriteField_u('f_addu',f_addu)
[1632]	605
[1806]	606	if (guide_zon) CALL guide_zonave_u(1,llm,f_addu)
	607	CALL guide_addfield_u(llm,f_addu,alpha_u)
[2028]	608	IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:),factt)
	609	IF (f_out) CALL guide_out("u",jjp1,llm,ucov(ijb_u:ije_u,:),factt)
[3995]	610	IF (f_out) THEN
	611	! Ehouarn: fill the gaps adequately...
	612	IF (ijbu>ijb_u) f_addu(ijb_u:ijbu-1,:)=0
	613	IF (ijeu<ije_u) f_addu(ijeu+1:ije_u,:)=0
	614	CALL guide_out("ucov",jjp1,llm,f_addu(ijb_u:ije_u,:)/factt,factt)
	615	ENDIF
[2036]	616	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	617	DO l=1,llm
	618	ucov(ijbu:ijeu,l)=ucov(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
	619	ENDDO
	620
[1632]	621	endif
	622
	623	if (guide_T) then
	624	if (guide_add) then
[2036]	625	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	626	DO l=1,llm
	627	f_addu(ijbu:ijeu,l)=(1.-tau)tgui1(ijbu:ijeu,l)+tautgui2(ijbu:ijeu,l)
	628	ENDDO
[1632]	629	else
[2036]	630	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	631	DO l=1,llm
	632	f_addu(ijbu:ijeu,l)=(1.-tau)tgui1(ijbu:ijeu,l)+tautgui2(ijbu:ijeu,l)-teta(ijbu:ijeu,l)
	633	ENDDO
[1632]	634	endif
[1806]	635	if (guide_zon) CALL guide_zonave_u(2,llm,f_addu)
	636	CALL guide_addfield_u(llm,f_addu,alpha_T)
[2028]	637	IF (f_out) CALL guide_out("teta",jjp1,llm,f_addu(:,:)/factt,factt)
[2036]	638	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	639	DO l=1,llm
	640	teta(ijbu:ijeu,l)=teta(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
	641	ENDDO
[1632]	642	endif
	643
	644	if (guide_P) then
	645	if (guide_add) then
[1806]	646	!$OMP MASTER
	647	f_addu(ijbu:ijeu,1)=(1.-tau)psgui1(ijbu:ijeu)+taupsgui2(ijbu:ijeu)
	648	!$OMP END MASTER
	649	!$OMP BARRIER
[1632]	650	else
[1806]	651	!$OMP MASTER
	652	f_addu(ijbu:ijeu,1)=(1.-tau)psgui1(ijbu:ijeu)+taupsgui2(ijbu:ijeu)-ps(ijbu:ijeu)
	653	!$OMP END MASTER
	654	!$OMP BARRIER
[1632]	655	endif
[1806]	656	if (guide_zon) CALL guide_zonave_u(2,1,f_addu(ijb_u:ije_u,1))
	657	CALL guide_addfield_u(1,f_addu(ijb_u:ije_u,1),alpha_P)
[2036]	658	! IF (f_out) CALL guide_out("ps",jjp1,1,f_addu(ijb_u:ije_u,1)/factt,factt)
[1806]	659	!$OMP MASTER
	660	ps(ijbu:ijeu)=ps(ijbu:ijeu)+f_addu(ijbu:ijeu,1)
	661	!$OMP END MASTER
	662	!$OMP BARRIER
[1632]	663	CALL pression_loc(ijnb_u,ap,bp,ps,p)
	664	CALL massdair_loc(p,masse)
[1806]	665	!$OMP BARRIER
[1632]	666	endif
	667
	668	if (guide_Q) then
	669	if (guide_add) then
[2036]	670	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	671	DO l=1,llm
	672	f_addu(ijbu:ijeu,l)=(1.-tau)qgui1(ijbu:ijeu,l)+tauqgui2(ijbu:ijeu,l)
	673	ENDDO
[1632]	674	else
[2036]	675	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	676	DO l=1,llm
	677	f_addu(ijbu:ijeu,l)=(1.-tau)qgui1(ijbu:ijeu,l)+tauqgui2(ijbu:ijeu,l)-q(ijbu:ijeu,l)
	678	ENDDO
[1632]	679	endif
[1806]	680	if (guide_zon) CALL guide_zonave_u(2,llm,f_addu)
	681	CALL guide_addfield_u(llm,f_addu,alpha_Q)
[2028]	682	IF (f_out) CALL guide_out("q",jjp1,llm,f_addu(:,:)/factt,factt)
[1806]	683
[2036]	684	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	685	DO l=1,llm
	686	q(ijbu:ijeu,l)=q(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
	687	ENDDO
[1632]	688	endif
	689
	690	if (guide_v) then
	691	if (guide_add) then
[2036]	692	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	693	DO l=1,llm
	694	f_addv(ijbv:ijev,l)=(1.-tau)vgui1(ijbv:ijev,l)+tauvgui2(ijbv:ijev,l)
	695	ENDDO
	696
[1632]	697	else
[2036]	698	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	699	DO l=1,llm
	700	f_addv(ijbv:ijev,l)=(1.-tau)vgui1(ijbv:ijev,l)+tauvgui2(ijbv:ijev,l)-vcov(ijbv:ijev,l)
	701	ENDDO
	702
[1632]	703	endif
[1806]	704
	705	if (guide_zon) CALL guide_zonave_v(2,jjm,llm,f_addv(ijb_v:ije_v,:))
[1632]	706
[1806]	707	CALL guide_addfield_v(llm,f_addv(ijb_v:ije_v,:),alpha_v)
[2028]	708	IF (f_out) CALL guide_out("v",jjm,llm,vcov(ijb_v:ije_v,:),factt)
	709	IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:),factt)
[3995]	710	IF (f_out) THEN
	711	! Ehouarn: Fill in the gaps adequately
	712	IF (ijbv>ijb_v) f_addv(ijb_v:ijbv-1,:)=0
	713	IF (ijev<ije_v) f_addv(ijev+1:ije_v,:)=0
	714	CALL guide_out("vcov",jjm,llm,f_addv(ijb_v:ije_v,:)/factt,factt)
	715	ENDIF
[1806]	716
[2036]	717	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	718	DO l=1,llm
	719	vcov(ijbv:ijev,l)=vcov(ijbv:ijev,l)+f_addv(ijbv:ijev,l)
	720	ENDDO
[1632]	721	endif
	722
	723	END SUBROUTINE guide_main
	724
	725
	726	SUBROUTINE guide_addfield_u(vsize,field,alpha)
	727	! field1=afield1+alphafield2
	728
	729	IMPLICIT NONE
	730	INCLUDE "dimensions.h"
	731	INCLUDE "paramet.h"
	732
	733	! input variables
	734	INTEGER, INTENT(IN) :: vsize
	735	REAL, DIMENSION(ijb_u:ije_u), INTENT(IN) :: alpha
	736	REAL, DIMENSION(ijb_u:ije_u,vsize), INTENT(INOUT) :: field
	737
	738	! Local variables
	739	INTEGER :: l
	740
[2036]	741	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	742	DO l=1,vsize
	743	field(ijbu:ijeu,l)=alpha(ijbu:ijeu)field(ijbu:ijeu,l)alpha_pcor(l)
	744	ENDDO
	745
	746	END SUBROUTINE guide_addfield_u
	747
	748
	749	SUBROUTINE guide_addfield_v(vsize,field,alpha)
	750	! field1=afield1+alphafield2
	751
	752	IMPLICIT NONE
	753	INCLUDE "dimensions.h"
	754	INCLUDE "paramet.h"
	755
	756	! input variables
	757	INTEGER, INTENT(IN) :: vsize
	758	REAL, DIMENSION(ijb_v:ije_v), INTENT(IN) :: alpha
	759	REAL, DIMENSION(ijb_v:ije_v,vsize), INTENT(INOUT) :: field
	760
	761	! Local variables
	762	INTEGER :: l
	763
[2036]	764	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	765	DO l=1,vsize
	766	field(ijbv:ijev,l)=alpha(ijbv:ijev)field(ijbv:ijev,l)alpha_pcor(l)
	767	ENDDO
	768
	769	END SUBROUTINE guide_addfield_v
	770
	771	!=======================================================================
	772
	773	SUBROUTINE guide_zonave_u(typ,vsize,field)
	774
[2597]	775	USE comconst_mod, ONLY: pi
	776
[1632]	777	IMPLICIT NONE
	778
	779	INCLUDE "dimensions.h"
	780	INCLUDE "paramet.h"
	781	INCLUDE "comgeom.h"
	782
	783	! input/output variables
	784	INTEGER, INTENT(IN) :: typ
	785	INTEGER, INTENT(IN) :: vsize
	786	REAL, DIMENSION(ijb_u:ije_u,vsize), INTENT(INOUT) :: field
	787
	788	! Local variables
	789	LOGICAL, SAVE :: first=.TRUE.
[1806]	790	!$OMP THREADPRIVATE(first)
	791
[1632]	792	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
[1806]	793	!$OMP THREADPRIVATE(imin,imax)
[1632]	794	INTEGER :: i,j,l,ij
	795	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
	796	REAL, DIMENSION (jjb_u:jje_u,vsize):: fieldm ! zon-averaged field
	797
	798	IF (first) THEN
	799	first=.FALSE.
	800	!Compute domain for averaging
	801	lond=rlonu*180./pi
	802	imin(1)=1;imax(1)=iip1;
	803	imin(2)=1;imax(2)=iip1;
	804	IF (guide_reg) THEN
	805	DO i=1,iim
[5084]	806	IF (lond(i).LT.lon_min_g) imin(1)=i
	807	IF (lond(i).LE.lon_max_g) imax(1)=i
[1632]	808	ENDDO
	809	lond=rlonv*180./pi
	810	DO i=1,iim
[5084]	811	IF (lond(i).LT.lon_min_g) imin(2)=i
	812	IF (lond(i).LE.lon_max_g) imax(2)=i
[1632]	813	ENDDO
	814	ENDIF
	815	ENDIF
	816
	817
[2036]	818	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	819	DO l=1,vsize
[1806]	820	fieldm(:,l)=0.
[1632]	821	! Compute zonal average
	822
	823	!correction bug ici
	824	! ---> a verifier
	825	! ym DO j=jjbv,jjev
	826	DO j=jjbu,jjeu
	827	DO i=imin(typ),imax(typ)
	828	ij=(j-1)*iip1+i
	829	fieldm(j,l)=fieldm(j,l)+field(ij,l)
	830	ENDDO
	831	ENDDO
	832	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
	833	! Compute forcing
	834	DO j=jjbu,jjeu
	835	DO i=1,iip1
	836	ij=(j-1)*iip1+i
	837	field(ij,l)=fieldm(j,l)
	838	ENDDO
	839	ENDDO
	840	ENDDO
	841
	842	END SUBROUTINE guide_zonave_u
	843
	844
	845	SUBROUTINE guide_zonave_v(typ,hsize,vsize,field)
	846
[2597]	847	USE comconst_mod, ONLY: pi
	848
[1632]	849	IMPLICIT NONE
	850
	851	INCLUDE "dimensions.h"
	852	INCLUDE "paramet.h"
	853	INCLUDE "comgeom.h"
	854
	855	! input/output variables
	856	INTEGER, INTENT(IN) :: typ
	857	INTEGER, INTENT(IN) :: vsize
	858	INTEGER, INTENT(IN) :: hsize
	859	REAL, DIMENSION(ijb_v:ije_v,vsize), INTENT(INOUT) :: field
	860
	861	! Local variables
	862	LOGICAL, SAVE :: first=.TRUE.
[1806]	863	!$OMP THREADPRIVATE(first)
[1632]	864	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
[1806]	865	!$OMP THREADPRIVATE(imin, imax)
[1632]	866	INTEGER :: i,j,l,ij
	867	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
	868	REAL, DIMENSION (jjb_v:jjev,vsize):: fieldm ! zon-averaged field
	869
	870	IF (first) THEN
	871	first=.FALSE.
	872	!Compute domain for averaging
	873	lond=rlonu*180./pi
	874	imin(1)=1;imax(1)=iip1;
	875	imin(2)=1;imax(2)=iip1;
	876	IF (guide_reg) THEN
	877	DO i=1,iim
[5084]	878	IF (lond(i).LT.lon_min_g) imin(1)=i
	879	IF (lond(i).LE.lon_max_g) imax(1)=i
[1632]	880	ENDDO
	881	lond=rlonv*180./pi
	882	DO i=1,iim
[5084]	883	IF (lond(i).LT.lon_min_g) imin(2)=i
	884	IF (lond(i).LE.lon_max_g) imax(2)=i
[1632]	885	ENDDO
	886	ENDIF
	887	ENDIF
	888
[2036]	889	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	890	DO l=1,vsize
	891	! Compute zonal average
[1806]	892	fieldm(:,l)=0.
[1632]	893	DO j=jjbv,jjev
	894	DO i=imin(typ),imax(typ)
	895	ij=(j-1)*iip1+i
	896	fieldm(j,l)=fieldm(j,l)+field(ij,l)
	897	ENDDO
	898	ENDDO
	899	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
	900	! Compute forcing
	901	DO j=jjbv,jjev
	902	DO i=1,iip1
	903	ij=(j-1)*iip1+i
	904	field(ij,l)=fieldm(j,l)
	905	ENDDO
	906	ENDDO
	907	ENDDO
	908
	909
	910	END SUBROUTINE guide_zonave_v
	911
	912	!=======================================================================
	913	SUBROUTINE guide_interp(psi,teta)
[2021]	914	use exner_hyb_loc_m, only: exner_hyb_loc
	915	use exner_milieu_loc_m, only: exner_milieu_loc
[1823]	916	USE parallel_lmdz
[1632]	917	USE mod_hallo
	918	USE Bands
[2597]	919	USE comconst_mod, ONLY: cpp, kappa
[2600]	920	USE comvert_mod, ONLY: preff, pressure_exner, bp, ap, disvert_type
[1632]	921	IMPLICIT NONE
	922
	923	include "dimensions.h"
	924	include "paramet.h"
	925	include "comgeom2.h"
	926
	927	REAL, DIMENSION (iip1,jjb_u:jje_u), INTENT(IN) :: psi ! Psol gcm
	928	REAL, DIMENSION (iip1,jjb_u:jje_u,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
	929
	930	LOGICAL, SAVE :: first=.TRUE.
[1806]	931	!$OMP THREADPRIVATE(first)
[1632]	932	! Variables pour niveaux pression:
[1806]	933	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: plnc1,plnc2 !niveaux pression guidage
	934	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: plunc,plsnc !niveaux pression modele
	935	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: plvnc !niveaux pression modele
	936	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: p ! pression intercouches
	937	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: pls, pext ! var intermediaire
	938	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: pbarx
	939	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: pbary
[1632]	940	! Variables pour fonction Exner (P milieu couche)
[2021]	941	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: pk
[1806]	942	REAL ,ALLOCATABLE, SAVE, DIMENSION (:,:) :: pks
[1632]	943	REAL :: unskap
	944	! Pression de vapeur saturante
[1806]	945	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:) :: qsat
[1632]	946	!Variables intermediaires interpolation
[1806]	947	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: zu1,zu2
	948	REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:) :: zv1,zv2
[1632]	949
	950	INTEGER :: i,j,l,ij
[3995]	951	CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
[1848]	952	TYPE(Request),SAVE :: Req
	953	!$OMP THREADPRIVATE(Req)
[3995]	954
	955	if (is_master) write(,)trim(modname)//': interpolate nudging variables'
[1632]	956	! -----------------------------------------------------------------
	957	! Calcul des niveaux de pression champs guidage (pour T et Q)
	958	! -----------------------------------------------------------------
[1806]	959	IF (first) THEN
	960	!$OMP MASTER
	961	ALLOCATE(plnc1(iip1,jjb_u:jje_u,nlevnc) )
	962	ALLOCATE(plnc2(iip1,jjb_u:jje_u,nlevnc) )
	963	ALLOCATE(plunc(iip1,jjb_u:jje_u,llm) )
	964	ALLOCATE(plsnc(iip1,jjb_u:jje_u,llm) )
	965	ALLOCATE(plvnc(iip1,jjb_v:jje_v,llm) )
	966	ALLOCATE(p(iip1,jjb_u:jje_u,llmp1) )
	967	ALLOCATE(pls(iip1,jjb_u:jje_u,llm) )
	968	ALLOCATE(pext(iip1,jjb_u:jje_u,llm) )
	969	ALLOCATE(pbarx(iip1,jjb_u:jje_u,llm) )
	970	ALLOCATE(pbary(iip1,jjb_v:jje_v,llm) )
	971	ALLOCATE(pk(iip1,jjb_u:jje_u,llm) )
	972	ALLOCATE(pks (iip1,jjb_u:jje_u) )
	973	ALLOCATE(qsat(ijb_u:ije_u,llm) )
	974	ALLOCATE(zu1(iip1,jjb_u:jje_u,llm) )
	975	ALLOCATE(zu2(iip1,jjb_u:jje_u,llm) )
	976	ALLOCATE(zv1(iip1,jjb_v:jje_v,llm) )
	977	ALLOCATE(zv2(iip1,jjb_v:jje_v,llm) )
	978	!$OMP END MASTER
	979	!$OMP BARRIER
	980	ENDIF
	981
	982
	983
	984
[5084]	985	IF (guide_plevs.EQ.0) THEN
[1806]	986	!$OMP DO
[1632]	987	DO l=1,nlevnc
	988	DO j=jjbu,jjeu
	989	DO i=1,iip1
	990	plnc2(i,j,l)=apnc(l)
	991	plnc1(i,j,l)=apnc(l)
	992	ENDDO
	993	ENDDO
	994	ENDDO
	995	ENDIF
	996
	997	if (first) then
	998	first=.FALSE.
[1806]	999	!$OMP MASTER
[3995]	1000	write(,)trim(modname)//' : check vertical level order'
	1001	write(,)trim(modname)//' LMDZ :'
[1632]	1002	do l=1,llm
[3995]	1003	write(,)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
[1632]	1004	+psi(1,jjeu)*(bp(l)+bp(l+1))/2.
	1005	enddo
[3995]	1006	write(,)trim(modname)//' nudging file :'
[1632]	1007	SELECT CASE (guide_plevs)
	1008	CASE (0)
	1009	do l=1,nlevnc
[3995]	1010	write(,)trim(modname)//' PL(',l,')=',plnc2(1,jjbu,l)
[1632]	1011	enddo
	1012	CASE (1)
	1013	DO l=1,nlevnc
[3995]	1014	write(,)trim(modname)//' PL(',l,')=',&
[4041]	1015	apnc(l)+bpnc(l)*psnat2(1,jjbu)
[3995]	1016	ENDDO
[1632]	1017	CASE (2)
	1018	do l=1,nlevnc
[3995]	1019	write(,)trim(modname)//' PL(',l,')=',pnat2(1,jjbu,l)
[1632]	1020	enddo
	1021	END SELECT
[3995]	1022	write(,)trim(modname)//' invert ordering: invert_p=',invert_p
[1632]	1023	if (guide_u) then
	1024	do l=1,nlevnc
[3995]	1025	write(,)trim(modname)//' U(',l,')=',unat2(1,jjbu,l)
[1632]	1026	enddo
	1027	endif
	1028	if (guide_T) then
	1029	do l=1,nlevnc
[3995]	1030	write(,)trim(modname)//' T(',l,')=',tnat2(1,jjbu,l)
[1632]	1031	enddo
	1032	endif
[1806]	1033	!$OMP END MASTER
[3995]	1034	endif ! of if (first)
[4042]	1035
[4054]	1036	if (guide_plevs /= 1 .or. guide_t .and. .not. guide_teta &
	1037	.or. guide_q .and. guide_hr) then
	1038	CALL pression_loc( ijnb_u, ap, bp, psi, p )
	1039	if (disvert_type==1) then
	1040	CALL exner_hyb_loc(ijnb_u,psi,p,pks,pk)
	1041	else ! we assume that we are in the disvert_type==2 case
	1042	CALL exner_milieu_loc(ijnb_u,psi,p,pks,pk)
	1043	endif
	1044	end if
	1045
[1632]	1046	! -----------------------------------------------------------------
	1047	! Calcul niveaux pression modele
	1048	! -----------------------------------------------------------------
	1049
	1050	! .... Calcul de pls , pression au milieu des couches ,en Pascals
[5084]	1051	IF (guide_plevs.EQ.1) THEN
[2036]	1052	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1053	DO l=1,llm
[2597]	1054	DO j=jjbu,jjeu
	1055	DO i =1, iip1
[1632]	1056	pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
[2028]	1057	ENDDO
	1058	ENDDO
[1632]	1059	ENDDO
	1060	ELSE
[2028]	1061	unskap=1./kappa
[1806]	1062	!$OMP BARRIER
[2036]	1063	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2028]	1064	DO l = 1, llm
	1065	DO j=jjbu,jjeu
[2597]	1066	DO i =1, iip1
	1067	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
	1068	ENDDO
[2028]	1069	ENDDO
	1070	ENDDO
[1632]	1071	ENDIF
	1072
	1073	! calcul des pressions pour les grilles u et v
[2036]	1074	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1075	do l=1,llm
	1076	do j=jjbu,jjeu
	1077	do i=1,iip1
	1078	pext(i,j,l)=pls(i,j,l)*aire(i,j)
	1079	enddo
	1080	enddo
	1081	enddo
	1082
	1083	CALL Register_Hallo_u(pext,llm,1,2,2,1,Req)
	1084	CALL SendRequest(Req)
[1806]	1085	!$OMP BARRIER
[1632]	1086	CALL WaitRequest(Req)
[1806]	1087	!$OMP BARRIER
[1632]	1088
	1089	call massbar_loc(pext, pbarx, pbary )
[1806]	1090	!$OMP BARRIER
[2036]	1091	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1092	do l=1,llm
	1093	do j=jjbu,jjeu
	1094	do i=1,iip1
	1095	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
	1096	plsnc(i,j,l)=pls(i,j,l)
	1097	enddo
	1098	enddo
	1099	enddo
[2036]	1100	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1101	do l=1,llm
	1102	do j=jjbv,jjev
	1103	do i=1,iip1
	1104	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
	1105	enddo
	1106	enddo
	1107	enddo
	1108
	1109	! -----------------------------------------------------------------
	1110	! Interpolation verticale champs guidage sur niveaux modele
	1111	! Conversion en variables gcm (ucov, vcov...)
	1112	! -----------------------------------------------------------------
	1113	if (guide_P) then
[1806]	1114	!$OMP MASTER
[1632]	1115	do j=jjbu,jjeu
	1116	do i=1,iim
	1117	ij=(j-1)*iip1+i
	1118	psgui1(ij)=psnat1(i,j)
	1119	psgui2(ij)=psnat2(i,j)
	1120	enddo
	1121	psgui1(iip1*j)=psnat1(1,j)
	1122	psgui2(iip1*j)=psnat2(1,j)
	1123	enddo
[1806]	1124	!$OMP END MASTER
	1125	!$OMP BARRIER
[1632]	1126	endif
	1127
	1128	IF (guide_T) THEN
	1129	! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]	1130	IF (guide_plevs.EQ.1) THEN
[1806]	1131	!$OMP DO
[1632]	1132	DO l=1,nlevnc
	1133	DO j=jjbu,jjeu
	1134	DO i=1,iip1
	1135	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
	1136	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
	1137	ENDDO
	1138	ENDDO
	1139	ENDDO
[5084]	1140	ELSE IF (guide_plevs.EQ.2) THEN
[1806]	1141	!$OMP DO
[1632]	1142	DO l=1,nlevnc
	1143	DO j=jjbu,jjeu
	1144	DO i=1,iip1
	1145	plnc2(i,j,l)=pnat2(i,j,l)
	1146	plnc1(i,j,l)=pnat1(i,j,l)
	1147	ENDDO
	1148	ENDDO
	1149	ENDDO
	1150	ENDIF
	1151
	1152	! Interpolation verticale
[1806]	1153	!$OMP MASTER
[1632]	1154	CALL pres2lev(tnat1(:,jjbu:jjeu,:),zu1(:,jjbu:jjeu,:),nlevnc,llm, &
	1155	plnc1(:,jjbu:jjeu,:),plsnc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
	1156	CALL pres2lev(tnat2(:,jjbu:jjeu,:),zu2(:,jjbu:jjeu,:),nlevnc,llm, &
	1157	plnc2(:,jjbu:jjeu,:),plsnc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
[1806]	1158	!$OMP END MASTER
	1159	!$OMP BARRIER
[1632]	1160	! Conversion en variables GCM
[2036]	1161	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1162	do l=1,llm
	1163	do j=jjbu,jjeu
	1164	IF (guide_teta) THEN
	1165	do i=1,iim
	1166	ij=(j-1)*iip1+i
	1167	tgui1(ij,l)=zu1(i,j,l)
	1168	tgui2(ij,l)=zu2(i,j,l)
	1169	enddo
	1170	ELSE
	1171	do i=1,iim
	1172	ij=(j-1)*iip1+i
	1173	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
	1174	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
	1175	enddo
	1176	ENDIF
	1177	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
	1178	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
	1179	enddo
	1180	if (pole_nord) then
	1181	do i=1,iip1
	1182	tgui1(i,l)=tgui1(1,l)
	1183	tgui2(i,l)=tgui2(1,l)
	1184	enddo
	1185	endif
	1186	if (pole_sud) then
	1187	do i=1,iip1
	1188	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
	1189	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
	1190	enddo
	1191	endif
	1192	enddo
	1193	ENDIF
	1194
	1195	IF (guide_Q) THEN
	1196	! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]	1197	IF (guide_plevs.EQ.1) THEN
[1806]	1198	!$OMP DO
[1632]	1199	DO l=1,nlevnc
	1200	DO j=jjbu,jjeu
	1201	DO i=1,iip1
	1202	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
	1203	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
	1204	ENDDO
	1205	ENDDO
	1206	ENDDO
[5084]	1207	ELSE IF (guide_plevs.EQ.2) THEN
[1806]	1208	!$OMP DO
[1632]	1209	DO l=1,nlevnc
	1210	DO j=jjbu,jjeu
	1211	DO i=1,iip1
	1212	plnc2(i,j,l)=pnat2(i,j,l)
	1213	plnc1(i,j,l)=pnat1(i,j,l)
	1214	ENDDO
	1215	ENDDO
	1216	ENDDO
	1217	ENDIF
	1218
	1219	! Interpolation verticale
[1806]	1220	!$OMP MASTER
[1632]	1221	CALL pres2lev(qnat1(:,jjbu:jjeu,:),zu1(:,jjbu:jjeu,:),nlevnc,llm, &
	1222	plnc1(:,jjbu:jjeu,:),plsnc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
	1223	CALL pres2lev(qnat2(:,jjbu:jjeu,:),zu2(:,jjbu:jjeu,:),nlevnc,llm, &
	1224	plnc2(:,jjbu:jjeu,:),plsnc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
[1806]	1225	!$OMP END MASTER
	1226	!$OMP BARRIER
[1632]	1227
	1228	! Conversion en variables GCM
	1229	! On suppose qu'on a la bonne variable dans le fichier de guidage:
	1230	! Hum.Rel si guide_hr, Hum.Spec. sinon.
[2036]	1231	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1232	do l=1,llm
	1233	do j=jjbu,jjeu
	1234	do i=1,iim
	1235	ij=(j-1)*iip1+i
	1236	qgui1(ij,l)=zu1(i,j,l)
	1237	qgui2(ij,l)=zu2(i,j,l)
	1238	enddo
	1239	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
	1240	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
	1241	enddo
	1242	if (pole_nord) then
	1243	do i=1,iip1
	1244	qgui1(i,l)=qgui1(1,l)
	1245	qgui2(i,l)=qgui2(1,l)
	1246	enddo
	1247	endif
[3547]	1248	if (pole_sud) then
[1632]	1249	do i=1,iip1
	1250	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
	1251	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
	1252	enddo
	1253	endif
	1254	enddo
	1255	IF (guide_hr) THEN
[2036]	1256	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]	1257	do l=1,llm
	1258	CALL q_sat(iip1jjnu,teta(:,jjbu:jjeu,l)pk(:,jjbu:jjeu,l)/cpp, &
	1259	plsnc(:,jjbu:jjeu,l),qsat(ijbu:ijeu,l))
	1260	qgui1(ijbu:ijeu,l)=qgui1(ijbu:ijeu,l)qsat(ijbu:ijeu,l)0.01 !hum. rel. en %
	1261	qgui2(ijbu:ijeu,l)=qgui2(ijbu:ijeu,l)qsat(ijbu:ijeu,l)0.01
	1262	enddo
	1263
[1632]	1264	ENDIF
	1265	ENDIF
	1266
	1267	IF (guide_u) THEN
	1268	! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]	1269	IF (guide_plevs.EQ.1) THEN
[1806]	1270	!$OMP DO
[1632]	1271	DO l=1,nlevnc
	1272	DO j=jjbu,jjeu
	1273	DO i=1,iim
	1274	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha1p2(i,j) &
	1275	& +psnat2(i+1,j)aire(i+1,j)alpha3p4(i+1,j))/aireu(i,j)
	1276	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha1p2(i,j) &
	1277	& +psnat1(i+1,j)aire(i+1,j)alpha3p4(i+1,j))/aireu(i,j)
	1278	ENDDO
	1279	plnc2(iip1,j,l)=plnc2(1,j,l)
	1280	plnc1(iip1,j,l)=plnc1(1,j,l)
	1281	ENDDO
	1282	ENDDO
[5084]	1283	ELSE IF (guide_plevs.EQ.2) THEN
[1806]	1284	!$OMP DO
[1632]	1285	DO l=1,nlevnc
	1286	DO j=jjbu,jjeu
	1287	DO i=1,iim
	1288	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha1p2(i,j) &
	1289	& +pnat2(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
	1290	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha1p2(i,j) &
	1291	& +pnat1(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
	1292	ENDDO
	1293	plnc2(iip1,j,l)=plnc2(1,j,l)
	1294	plnc1(iip1,j,l)=plnc1(1,j,l)
	1295	ENDDO
	1296	ENDDO
	1297	ENDIF
	1298
	1299	! Interpolation verticale
[1806]	1300	!$OMP MASTER
[1632]	1301	CALL pres2lev(unat1(:,jjbu:jjeu,:),zu1(:,jjbu:jjeu,:),nlevnc,llm, &
	1302	plnc1(:,jjbu:jjeu,:),plunc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
	1303	CALL pres2lev(unat2(:,jjbu:jjeu,:),zu2(:,jjbu:jjeu,:),nlevnc,llm, &
	1304	plnc2(:,jjbu:jjeu,:),plunc(:,jjbu:jjeu,:),iip1,jjnu,invert_p)
[1806]	1305	!$OMP END MASTER
	1306	!$OMP BARRIER
[1632]	1307
	1308	! Conversion en variables GCM
[2036]	1309	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1310	do l=1,llm
	1311	do j=jjbu,jjeu
	1312	do i=1,iim
	1313	ij=(j-1)*iip1+i
	1314	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
	1315	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
	1316	enddo
	1317	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
	1318	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
	1319	enddo
	1320	if (pole_nord) then
	1321	do i=1,iip1
	1322	ugui1(i,l)=0.
	1323	ugui2(i,l)=0.
	1324	enddo
	1325	endif
	1326	if (pole_sud) then
	1327	do i=1,iip1
	1328	ugui1(ip1jm+i,l)=0.
	1329	ugui2(ip1jm+i,l)=0.
	1330	enddo
	1331	endif
	1332	enddo
	1333	ENDIF
	1334
	1335	IF (guide_v) THEN
	1336	! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]	1337	IF (guide_plevs.EQ.1) THEN
[2597]	1338	CALL Register_Hallo_u(psnat1,1,1,2,2,1,Req)
	1339	CALL Register_Hallo_u(psnat2,1,1,2,2,1,Req)
	1340	CALL SendRequest(Req)
[1806]	1341	!$OMP BARRIER
[2597]	1342	CALL WaitRequest(Req)
[1806]	1343	!$OMP BARRIER
	1344	!$OMP DO
[1632]	1345	DO l=1,nlevnc
	1346	DO j=jjbv,jjev
	1347	DO i=1,iip1
	1348	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha2p3(i,j) &
	1349	& +psnat2(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
	1350	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha2p3(i,j) &
	1351	& +psnat1(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
	1352	ENDDO
	1353	ENDDO
	1354	ENDDO
[5084]	1355	ELSE IF (guide_plevs.EQ.2) THEN
[2597]	1356	CALL Register_Hallo_u(pnat1,llm,1,2,2,1,Req)
	1357	CALL Register_Hallo_u(pnat2,llm,1,2,2,1,Req)
	1358	CALL SendRequest(Req)
[1806]	1359	!$OMP BARRIER
[2597]	1360	CALL WaitRequest(Req)
[1806]	1361	!$OMP BARRIER
	1362	!$OMP DO
[1632]	1363	DO l=1,nlevnc
	1364	DO j=jjbv,jjev
	1365	DO i=1,iip1
	1366	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha2p3(i,j) &
	1367	& +pnat2(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
	1368	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha2p3(i,j) &
	1369	& +pnat1(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
	1370	ENDDO
	1371	ENDDO
	1372	ENDDO
	1373	ENDIF
	1374	! Interpolation verticale
[1806]	1375
	1376	!$OMP MASTER
[1632]	1377	CALL pres2lev(vnat1(:,jjbv:jjev,:),zv1(:,jjbv:jjev,:),nlevnc,llm, &
	1378	plnc1(:,jjbv:jjev,:),plvnc(:,jjbv:jjev,:),iip1,jjnv,invert_p)
	1379	CALL pres2lev(vnat2(:,jjbv:jjev,:),zv2(:,jjbv:jjev,:),nlevnc,llm, &
	1380	plnc2(:,jjbv:jjev,:),plvnc(:,jjbv:jjev,:),iip1,jjnv,invert_p)
[1806]	1381	!$OMP END MASTER
	1382	!$OMP BARRIER
[1632]	1383	! Conversion en variables GCM
[2036]	1384	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1385	do l=1,llm
	1386	do j=jjbv,jjev
	1387	do i=1,iim
	1388	ij=(j-1)*iip1+i
	1389	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
	1390	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
	1391	enddo
	1392	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
	1393	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
	1394	enddo
	1395	enddo
	1396	ENDIF
	1397
	1398
	1399	END SUBROUTINE guide_interp
	1400
	1401	!=======================================================================
[1806]	1402	SUBROUTINE tau2alpha(typ,pim,jjb,jje,factt,taumin,taumax,alpha)
[1632]	1403
	1404	! Calcul des constantes de rappel alpha (=1/tau)
	1405
[2597]	1406	use comconst_mod, only: pi
[2598]	1407	use serre_mod, only: clat, clon, grossismx, grossismy
[2597]	1408
[1632]	1409	implicit none
	1410
	1411	include "dimensions.h"
	1412	include "paramet.h"
	1413	include "comgeom2.h"
	1414
	1415	! input arguments :
	1416	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
[1806]	1417	INTEGER, INTENT(IN) :: pim ! dimensions en lon
	1418	INTEGER, INTENT(IN) :: jjb,jje ! dimensions en lat
[1632]	1419	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
	1420	REAL, INTENT(IN) :: taumin,taumax
	1421	! output arguments:
[1806]	1422	REAL, DIMENSION(pim,jjb:jje), INTENT(OUT) :: alpha
[1632]	1423
	1424	! local variables:
	1425	LOGICAL, SAVE :: first=.TRUE.
	1426	REAL, SAVE :: gamma,dxdy_min,dxdy_max
	1427	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
	1428	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
	1429	REAL, DIMENSION (iip1,jjm) :: dxdyv
	1430	real dxdy_
	1431	real zlat,zlon
	1432	real alphamin,alphamax,xi
	1433	integer i,j,ilon,ilat
[3995]	1434	character(len=20),parameter :: modname="tau2alpha"
[1632]	1435
	1436
	1437	alphamin=factt/taumax
	1438	alphamax=factt/taumin
	1439	IF (guide_reg.OR.guide_add) THEN
	1440	alpha=alphamax
	1441	!-----------------------------------------------------------------------
	1442	! guide_reg: alpha=alpha_min dans region, 0. sinon.
	1443	!-----------------------------------------------------------------------
	1444	IF (guide_reg) THEN
[1806]	1445	do j=jjb,jje
[1632]	1446	do i=1,pim
[5084]	1447	if (typ.eq.2) then
[1632]	1448	zlat=rlatu(j)*180./pi
	1449	zlon=rlonu(i)*180./pi
[5084]	1450	elseif (typ.eq.1) then
[1632]	1451	zlat=rlatu(j)*180./pi
	1452	zlon=rlonv(i)*180./pi
[5084]	1453	elseif (typ.eq.3) then
[1632]	1454	zlat=rlatv(j)*180./pi
	1455	zlon=rlonv(i)*180./pi
	1456	endif
	1457	alpha(i,j)=alphamax/16.* &
	1458	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
	1459	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
	1460	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
	1461	(1.+tanh((lon_max_g-zlon)/tau_lon))
	1462	enddo
	1463	enddo
	1464	ENDIF
	1465	ELSE
	1466	!-----------------------------------------------------------------------
	1467	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
	1468	!-----------------------------------------------------------------------
	1469	!Calcul de l'aire des mailles
	1470	do j=2,jjm
	1471	do i=2,iip1
	1472	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
	1473	enddo
	1474	zdx(1,j)=zdx(iip1,j)
	1475	enddo
	1476	do j=2,jjm
	1477	do i=1,iip1
	1478	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
	1479	enddo
	1480	enddo
	1481	do i=1,iip1
	1482	zdx(i,1)=zdx(i,2)
	1483	zdx(i,jjp1)=zdx(i,jjm)
	1484	zdy(i,1)=zdy(i,2)
	1485	zdy(i,jjp1)=zdy(i,jjm)
	1486	enddo
	1487	do j=1,jjp1
	1488	do i=1,iip1
	1489	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
	1490	enddo
	1491	enddo
[5084]	1492	IF (typ.EQ.2) THEN
[1632]	1493	do j=1,jjp1
	1494	do i=1,iim
	1495	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
	1496	enddo
	1497	dxdyu(iip1,j)=dxdyu(1,j)
	1498	enddo
	1499	ENDIF
[5084]	1500	IF (typ.EQ.3) THEN
[1632]	1501	do j=1,jjm
	1502	do i=1,iip1
	1503	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
	1504	enddo
	1505	enddo
	1506	ENDIF
	1507	! Premier appel: calcul des aires min et max et de gamma.
	1508	IF (first) THEN
	1509	first=.FALSE.
	1510	! coordonnees du centre du zoom
	1511	CALL coordij(clon,clat,ilon,ilat)
	1512	! aire de la maille au centre du zoom
	1513	dxdy_min=dxdys(ilon,ilat)
	1514	! dxdy maximale de la maille
	1515	dxdy_max=0.
	1516	do j=1,jjp1
	1517	do i=1,iip1
	1518	dxdy_max=max(dxdy_max,dxdys(i,j))
	1519	enddo
	1520	enddo
	1521	! Calcul de gamma
[5084]	1522	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[3995]	1523	write(,)trim(modname)//' ATTENTION modele peu zoome'
	1524	write(,)trim(modname)//' ATTENTION on prend une constante de guidage cste'
	1525	gamma=0.
[1632]	1526	else
[3995]	1527	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
	1528	write(,)trim(modname)//' gamma=',gamma
[5084]	1529	if (gamma.lt.1.e-5) then
[3995]	1530	write(,)trim(modname)//' gamma =',gamma,'<1e-5'
[4469]	1531	CALL abort_gcm("guide_loc_mod","stopped",1)
[3995]	1532	endif
	1533	gamma=log(0.5)/log(gamma)
	1534	if (gamma4) then
	1535	gamma=min(gamma,4.)
	1536	endif
	1537	write(,)trim(modname)//' gamma=',gamma
[1632]	1538	endif
	1539	ENDIF !first
	1540
[1806]	1541	do j=jjb,jje
[1632]	1542	do i=1,pim
[5084]	1543	if (typ.eq.1) then
[1632]	1544	dxdy_=dxdys(i,j)
	1545	zlat=rlatu(j)*180./pi
[5084]	1546	elseif (typ.eq.2) then
[1632]	1547	dxdy_=dxdyu(i,j)
	1548	zlat=rlatu(j)*180./pi
[5084]	1549	elseif (typ.eq.3) then
[1632]	1550	dxdy_=dxdyv(i,j)
	1551	zlat=rlatv(j)*180./pi
	1552	endif
[5084]	1553	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[1632]	1554	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
	1555	alpha(i,j)=alphamin
	1556	else
	1557	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
	1558	xi=min(xi,1.)
[5084]	1559	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
[1632]	1560	alpha(i,j)=xialphamin+(1.-xi)alphamax
	1561	else
	1562	alpha(i,j)=0.
	1563	endif
	1564	endif
	1565	enddo
	1566	enddo
	1567	ENDIF ! guide_reg
	1568
[2134]	1569	if (.not. guide_add) alpha = 1. - exp(- alpha)
	1570
[1632]	1571	END SUBROUTINE tau2alpha
	1572
	1573	!=======================================================================
	1574	SUBROUTINE guide_read(timestep)
	1575
	1576	IMPLICIT NONE
	1577
[5084]	1578	include "netcdf.inc"
[3995]	1579	include "dimensions.h"
	1580	include "paramet.h"
[1632]	1581
	1582	INTEGER, INTENT(IN) :: timestep
	1583
	1584	LOGICAL, SAVE :: first=.TRUE.
	1585	! Identification fichiers et variables NetCDF:
	1586	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
	1587	INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
[3984]	1588	INTEGER :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
[1632]	1589	! Variables auxiliaires NetCDF:
	1590	INTEGER, DIMENSION(4) :: start,count
	1591	INTEGER :: status,rcode
	1592	CHARACTER (len = 80) :: abort_message
	1593	CHARACTER (len = 20) :: modname = 'guide_read'
[3984]	1594	CHARACTER (len = 20) :: namedim
[1632]	1595	abort_message='pb in guide_read'
	1596
	1597	! -----------------------------------------------------------------
	1598	! Premier appel: initialisation de la lecture des fichiers
	1599	! -----------------------------------------------------------------
	1600	if (first) then
	1601	ncidpl=-99
[3995]	1602	write(,),trim(modname)//': opening nudging files '
[1632]	1603	! Ap et Bp si Niveaux de pression hybrides
[5084]	1604	if (guide_plevs.EQ.1) then
[3995]	1605	write(,),trim(modname)//' Reading nudging on model levels'
[1632]	1606	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5084]	1607	IF (rcode.NE.NF_NOERR) THEN
[3995]	1608	abort_message='Nudging: error -> no file apbp.nc'
[1902]	1609	CALL abort_gcm(modname,abort_message,1)
	1610	ENDIF
[1632]	1611	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[5084]	1612	IF (rcode.NE.NF_NOERR) THEN
[3995]	1613	abort_message='Nudging: error -> no AP variable in file apbp.nc'
[1902]	1614	CALL abort_gcm(modname,abort_message,1)
	1615	ENDIF
[1632]	1616	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[5084]	1617	IF (rcode.NE.NF_NOERR) THEN
[3995]	1618	abort_message='Nudging: error -> no BP variable in file apbp.nc'
[1902]	1619	CALL abort_gcm(modname,abort_message,1)
	1620	ENDIF
[3995]	1621	write(,),trim(modname)//' ncidpl,varidap',ncidpl,varidap
[1632]	1622	endif
[3984]	1623
[1632]	1624	! Pression si guidage sur niveaux P variables
[5084]	1625	if (guide_plevs.EQ.2) then
[1632]	1626	rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5084]	1627	IF (rcode.NE.NF_NOERR) THEN
[3995]	1628	abort_message='Nudging: error -> no file P.nc'
[1902]	1629	CALL abort_gcm(modname,abort_message,1)
	1630	ENDIF
[1632]	1631	rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[5084]	1632	IF (rcode.NE.NF_NOERR) THEN
[3995]	1633	abort_message='Nudging: error -> no PRES variable in file P.nc'
[1902]	1634	CALL abort_gcm(modname,abort_message,1)
	1635	ENDIF
[3995]	1636	write(,),trim(modname)//' ncidp,varidp',ncidp,varidp
[5084]	1637	if (ncidpl.eq.-99) ncidpl=ncidp
[1632]	1638	endif
[3984]	1639
[1632]	1640	! Vent zonal
	1641	if (guide_u) then
	1642	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5084]	1643	IF (rcode.NE.NF_NOERR) THEN
[3995]	1644	abort_message='Nudging: error -> no file u.nc'
[1902]	1645	CALL abort_gcm(modname,abort_message,1)
	1646	ENDIF
[1632]	1647	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[5084]	1648	IF (rcode.NE.NF_NOERR) THEN
[3995]	1649	abort_message='Nudging: error -> no UWND variable in file u.nc'
[1902]	1650	CALL abort_gcm(modname,abort_message,1)
	1651	ENDIF
[3995]	1652	write(,),trim(modname)//' ncidu,varidu',ncidu,varidu
[5084]	1653	if (ncidpl.eq.-99) ncidpl=ncidu
[3984]	1654
	1655
	1656	status=NF90_INQ_DIMID(ncidu, "LONU", dimid)
	1657	status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
[5084]	1658	IF (lendim .NE. iip1) THEN
[3995]	1659	abort_message='dimension LONU different from iip1 in u.nc'
[3984]	1660	CALL abort_gcm(modname,abort_message,1)
	1661	ENDIF
	1662
	1663	status=NF90_INQ_DIMID(ncidu, "LATU", dimid)
	1664	status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
[5084]	1665	IF (lendim .NE. jjp1) THEN
[3995]	1666	abort_message='dimension LATU different from jjp1 in u.nc'
[3984]	1667	CALL abort_gcm(modname,abort_message,1)
	1668	ENDIF
	1669
[1632]	1670	endif
[3984]	1671
[1632]	1672	! Vent meridien
	1673	if (guide_v) then
	1674	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5084]	1675	IF (rcode.NE.NF_NOERR) THEN
[3995]	1676	abort_message='Nudging: error -> no file v.nc'
[1902]	1677	CALL abort_gcm(modname,abort_message,1)
	1678	ENDIF
[1632]	1679	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[5084]	1680	IF (rcode.NE.NF_NOERR) THEN
[3995]	1681	abort_message='Nudging: error -> no VWND variable in file v.nc'
[1902]	1682	CALL abort_gcm(modname,abort_message,1)
	1683	ENDIF
[3995]	1684	write(,),trim(modname)//' ncidv,varidv',ncidv,varidv
[5084]	1685	if (ncidpl.eq.-99) ncidpl=ncidv
[3984]	1686
	1687	status=NF90_INQ_DIMID(ncidv, "LONV", dimid)
	1688	status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
	1689
[5084]	1690	IF (lendim .NE. iip1) THEN
[3995]	1691	abort_message='dimension LONV different from iip1 in v.nc'
[3984]	1692	CALL abort_gcm(modname,abort_message,1)
	1693	ENDIF
	1694
	1695
	1696	status=NF90_INQ_DIMID(ncidv, "LATV", dimid)
	1697	status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
[5084]	1698	IF (lendim .NE. jjm) THEN
[3995]	1699	abort_message='dimension LATV different from jjm in v.nc'
[3984]	1700	CALL abort_gcm(modname,abort_message,1)
	1701	ENDIF
	1702
	1703	endif
	1704
[1632]	1705	! Temperature
	1706	if (guide_T) then
	1707	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5084]	1708	IF (rcode.NE.NF_NOERR) THEN
[3995]	1709	abort_message='Nudging: error -> no file T.nc'
[1902]	1710	CALL abort_gcm(modname,abort_message,1)
	1711	ENDIF
[1632]	1712	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[5084]	1713	IF (rcode.NE.NF_NOERR) THEN
[3995]	1714	abort_message='Nudging: error -> no AIR variable in file T.nc'
[1902]	1715	CALL abort_gcm(modname,abort_message,1)
	1716	ENDIF
[3995]	1717	write(,),trim(modname)//' ncidT,varidT',ncidt,varidt
[5084]	1718	if (ncidpl.eq.-99) ncidpl=ncidt
[3984]	1719
	1720	status=NF90_INQ_DIMID(ncidt, "LONV", dimid)
	1721	status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
[5084]	1722	IF (lendim .NE. iip1) THEN
[3995]	1723	abort_message='dimension LONV different from iip1 in T.nc'
[3984]	1724	CALL abort_gcm(modname,abort_message,1)
	1725	ENDIF
	1726
	1727	status=NF90_INQ_DIMID(ncidt, "LATU", dimid)
	1728	status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
[5084]	1729	IF (lendim .NE. jjp1) THEN
[3995]	1730	abort_message='dimension LATU different from jjp1 in T.nc'
[3984]	1731	CALL abort_gcm(modname,abort_message,1)
	1732	ENDIF
	1733
[1632]	1734	endif
[3984]	1735
[1632]	1736	! Humidite
	1737	if (guide_Q) then
	1738	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5084]	1739	IF (rcode.NE.NF_NOERR) THEN
[3995]	1740	abort_message='Nudging: error -> no file hur.nc'
[1902]	1741	CALL abort_gcm(modname,abort_message,1)
	1742	ENDIF
[1632]	1743	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[5084]	1744	IF (rcode.NE.NF_NOERR) THEN
[3995]	1745	abort_message='Nudging: error -> no RH variable in file hur.nc'
[1902]	1746	CALL abort_gcm(modname,abort_message,1)
	1747	ENDIF
[3995]	1748	write(,),trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
[5084]	1749	if (ncidpl.eq.-99) ncidpl=ncidQ
[3984]	1750
	1751
	1752	status=NF90_INQ_DIMID(ncidQ, "LONV", dimid)
	1753	status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
[5084]	1754	IF (lendim .NE. iip1) THEN
[3995]	1755	abort_message='dimension LONV different from iip1 in hur.nc'
[3984]	1756	CALL abort_gcm(modname,abort_message,1)
	1757	ENDIF
	1758
	1759	status=NF90_INQ_DIMID(ncidQ, "LATU", dimid)
	1760	status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
[5084]	1761	IF (lendim .NE. jjp1) THEN
[3995]	1762	abort_message='dimension LATU different from jjp1 in hur.nc'
[3984]	1763	CALL abort_gcm(modname,abort_message,1)
	1764	ENDIF
	1765
	1766
[1632]	1767	endif
	1768	! Pression de surface
[5084]	1769	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1632]	1770	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5084]	1771	IF (rcode.NE.NF_NOERR) THEN
[3995]	1772	abort_message='Nudging: error -> no file ps.nc'
[1902]	1773	CALL abort_gcm(modname,abort_message,1)
	1774	ENDIF
[1632]	1775	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[5084]	1776	IF (rcode.NE.NF_NOERR) THEN
[3995]	1777	abort_message='Nudging: error -> no SP variable in file ps.nc'
[1902]	1778	CALL abort_gcm(modname,abort_message,1)
	1779	ENDIF
[3995]	1780	write(,),trim(modname)//' ncidps,varidps',ncidps,varidps
[1632]	1781	endif
	1782	! Coordonnee verticale
[5084]	1783	if (guide_plevs.EQ.0) then
[1632]	1784	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]	1785	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[3995]	1786	write(,),trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1632]	1787	endif
	1788	! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]	1789	IF (guide_plevs.EQ.1) THEN
	1790	#ifdef NC_DOUBLE
	1791	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	1792	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	1793	#else
	1794	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	1795	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	1796	#endif
	1797	ELSEIF (guide_plevs.EQ.0) THEN
	1798	#ifdef NC_DOUBLE
	1799	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	1800	#else
	1801	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	1802	#endif
[3795]	1803	!FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
	1804	IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
[1632]	1805	bpnc(:)=0.
	1806	ENDIF
	1807	first=.FALSE.
	1808	ENDIF ! (first)
	1809
	1810	! -----------------------------------------------------------------
	1811	! lecture des champs u, v, T, Q, ps
	1812	! -----------------------------------------------------------------
	1813
	1814	! dimensions pour les champs scalaires et le vent zonal
	1815	start(1)=1
	1816	start(2)=jjb_u
	1817	start(3)=1
	1818	start(4)=timestep
	1819
	1820	count(1)=iip1
	1821	count(2)=jjnb_u
	1822	count(3)=nlevnc
	1823	count(4)=1
	1824
	1825	IF (invert_y) start(2)=jjp1-jje_u+1
	1826	! Pression
[5084]	1827	if (guide_plevs.EQ.2) then
	1828	#ifdef NC_DOUBLE
	1829	status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,pnat2)
	1830	#else
	1831	status=NF_GET_VARA_REAL(ncidp,varidp,start,count,pnat2)
	1832	#endif
[1632]	1833	IF (invert_y) THEN
	1834	! PRINT*,"Invertion impossible actuellement"
	1835	! CALL abort_gcm(modname,abort_message,1)
	1836	CALL invert_lat(iip1,jjnb_u,nlevnc,pnat2)
	1837	ENDIF
	1838	endif
	1839
	1840	! Vent zonal
	1841	if (guide_u) then
[5084]	1842	#ifdef NC_DOUBLE
	1843	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
	1844	#else
	1845	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
	1846	#endif
[1632]	1847	IF (invert_y) THEN
	1848	! PRINT*,"Invertion impossible actuellement"
	1849	! CALL abort_gcm(modname,abort_message,1)
	1850	CALL invert_lat(iip1,jjnb_u,nlevnc,unat2)
	1851	ENDIF
	1852
	1853	endif
	1854
[2028]	1855
[1632]	1856	! Temperature
	1857	if (guide_T) then
[5084]	1858	#ifdef NC_DOUBLE
	1859	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
	1860	#else
	1861	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
	1862	#endif
[1632]	1863	IF (invert_y) THEN
	1864	! PRINT*,"Invertion impossible actuellement"
	1865	! CALL abort_gcm(modname,abort_message,1)
	1866	CALL invert_lat(iip1,jjnb_u,nlevnc,tnat2)
	1867	ENDIF
	1868	endif
	1869
	1870	! Humidite
	1871	if (guide_Q) then
[5084]	1872	#ifdef NC_DOUBLE
	1873	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
	1874	#else
	1875	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
	1876	#endif
[1632]	1877	IF (invert_y) THEN
	1878	! PRINT*,"Invertion impossible actuellement"
	1879	! CALL abort_gcm(modname,abort_message,1)
	1880	CALL invert_lat(iip1,jjnb_u,nlevnc,qnat2)
	1881	ENDIF
	1882
	1883	endif
	1884
	1885	! Vent meridien
	1886	if (guide_v) then
	1887	start(2)=jjb_v
	1888	count(2)=jjnb_v
	1889	IF (invert_y) start(2)=jjm-jje_v+1
	1890
[5084]	1891	#ifdef NC_DOUBLE
	1892	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
	1893	#else
	1894	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
	1895	#endif
[1632]	1896	IF (invert_y) THEN
	1897	! PRINT*,"Invertion impossible actuellement"
	1898	! CALL abort_gcm(modname,abort_message,1)
	1899	CALL invert_lat(iip1,jjnb_v,nlevnc,vnat2)
	1900	ENDIF
	1901	endif
	1902
	1903	! Pression de surface
[5084]	1904	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1632]	1905	start(2)=jjb_u
	1906	start(3)=timestep
	1907	start(4)=0
	1908	count(2)=jjnb_u
	1909	count(3)=1
	1910	count(4)=0
	1911	IF (invert_y) start(2)=jjp1-jje_u+1
[5084]	1912	#ifdef NC_DOUBLE
	1913	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
	1914	#else
	1915	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
	1916	#endif
[1632]	1917	IF (invert_y) THEN
	1918	! PRINT*,"Invertion impossible actuellement"
	1919	! CALL abort_gcm(modname,abort_message,1)
	1920	CALL invert_lat(iip1,jjnb_u,1,psnat2)
	1921	ENDIF
	1922	endif
	1923
	1924	END SUBROUTINE guide_read
	1925
	1926	!=======================================================================
	1927	SUBROUTINE guide_read2D(timestep)
	1928
	1929	IMPLICIT NONE
	1930
[5084]	1931	include "netcdf.inc"
[3995]	1932	include "dimensions.h"
	1933	include "paramet.h"
[1632]	1934
	1935	INTEGER, INTENT(IN) :: timestep
	1936
	1937	LOGICAL, SAVE :: first=.TRUE.
	1938	! Identification fichiers et variables NetCDF:
	1939	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
	1940	INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
	1941	INTEGER :: ncidpl,varidpl,varidap,varidbp
	1942	! Variables auxiliaires NetCDF:
	1943	INTEGER, DIMENSION(4) :: start,count
	1944	INTEGER :: status,rcode
	1945	! Variables for 3D extension:
	1946	REAL, DIMENSION (jjb_u:jje_u,llm) :: zu
	1947	REAL, DIMENSION (jjb_v:jje_v,llm) :: zv
	1948	INTEGER :: i
	1949	CHARACTER (len = 80) :: abort_message
[1906]	1950	CHARACTER (len = 20) :: modname = 'guide_read2D'
[1632]	1951	abort_message='pb in guide_read2D'
	1952
	1953	! -----------------------------------------------------------------
	1954	! Premier appel: initialisation de la lecture des fichiers
	1955	! -----------------------------------------------------------------
	1956	if (first) then
	1957	ncidpl=-99
[3995]	1958	write(,)trim(modname)//' : opening nudging files '
[1632]	1959	! Ap et Bp si niveaux de pression hybrides
[5084]	1960	if (guide_plevs.EQ.1) then
[3995]	1961	write(,)trim(modname)//' Reading nudging on model levels'
	1962	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5084]	1963	IF (rcode.NE.NF_NOERR) THEN
[3995]	1964	abort_message='Nudging: error -> no file apbp.nc'
	1965	CALL abort_gcm(modname,abort_message,1)
	1966	ENDIF
	1967	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[5084]	1968	IF (rcode.NE.NF_NOERR) THEN
[3995]	1969	abort_message='Nudging: error -> no AP variable in file apbp.nc'
	1970	CALL abort_gcm(modname,abort_message,1)
	1971	ENDIF
	1972	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[5084]	1973	IF (rcode.NE.NF_NOERR) THEN
[3995]	1974	abort_message='Nudging: error -> no BP variable in file apbp.nc'
	1975	CALL abort_gcm(modname,abort_message,1)
	1976	ENDIF
	1977	write(,)trim(modname)//'ncidpl,varidap',ncidpl,varidap
[1632]	1978	endif
	1979	! Pression
[5084]	1980	if (guide_plevs.EQ.2) then
[3995]	1981	rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5084]	1982	IF (rcode.NE.NF_NOERR) THEN
[3995]	1983	abort_message='Nudging: error -> no file P.nc'
	1984	CALL abort_gcm(modname,abort_message,1)
	1985	ENDIF
	1986	rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[5084]	1987	IF (rcode.NE.NF_NOERR) THEN
[3995]	1988	abort_message='Nudging: error -> no PRES variable in file P.nc'
	1989	CALL abort_gcm(modname,abort_message,1)
	1990	ENDIF
	1991	write(,)trim(modname)//' ncidp,varidp',ncidp,varidp
[5084]	1992	if (ncidpl.eq.-99) ncidpl=ncidp
[1632]	1993	endif
	1994	! Vent zonal
	1995	if (guide_u) then
[3995]	1996	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5084]	1997	IF (rcode.NE.NF_NOERR) THEN
[3995]	1998	abort_message='Nudging: error -> no file u.nc'
	1999	CALL abort_gcm(modname,abort_message,1)
	2000	ENDIF
	2001	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[5084]	2002	IF (rcode.NE.NF_NOERR) THEN
[3995]	2003	abort_message='Nudging: error -> no UWND variable in file u.nc'
	2004	CALL abort_gcm(modname,abort_message,1)
	2005	ENDIF
	2006	write(,)trim(modname)//' ncidu,varidu',ncidu,varidu
[5084]	2007	if (ncidpl.eq.-99) ncidpl=ncidu
[1632]	2008	endif
[2028]	2009
[1632]	2010	! Vent meridien
	2011	if (guide_v) then
[3995]	2012	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5084]	2013	IF (rcode.NE.NF_NOERR) THEN
[3995]	2014	abort_message='Nudging: error -> no file v.nc'
	2015	CALL abort_gcm(modname,abort_message,1)
	2016	ENDIF
	2017	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[5084]	2018	IF (rcode.NE.NF_NOERR) THEN
[3995]	2019	abort_message='Nudging: error -> no VWND variable in file v.nc'
	2020	CALL abort_gcm(modname,abort_message,1)
	2021	ENDIF
	2022	write(,)trim(modname)//' ncidv,varidv',ncidv,varidv
[5084]	2023	if (ncidpl.eq.-99) ncidpl=ncidv
[3995]	2024	endif
[1632]	2025	! Temperature
	2026	if (guide_T) then
[3995]	2027	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5084]	2028	IF (rcode.NE.NF_NOERR) THEN
[3995]	2029	abort_message='Nudging: error -> no file T.nc'
	2030	CALL abort_gcm(modname,abort_message,1)
	2031	ENDIF
	2032	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[5084]	2033	IF (rcode.NE.NF_NOERR) THEN
[3995]	2034	abort_message='Nudging: error -> no AIR variable in file T.nc'
	2035	CALL abort_gcm(modname,abort_message,1)
	2036	ENDIF
	2037	write(,)trim(modname)//' ncidT,varidT',ncidt,varidt
[5084]	2038	if (ncidpl.eq.-99) ncidpl=ncidt
[1632]	2039	endif
	2040	! Humidite
	2041	if (guide_Q) then
[3995]	2042	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5084]	2043	IF (rcode.NE.NF_NOERR) THEN
[3995]	2044	abort_message='Nudging: error -> no file hur.nc'
	2045	CALL abort_gcm(modname,abort_message,1)
	2046	ENDIF
	2047	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[5084]	2048	IF (rcode.NE.NF_NOERR) THEN
[3995]	2049	abort_message='Nudging: error -> no RH,variable in file hur.nc'
	2050	CALL abort_gcm(modname,abort_message,1)
	2051	ENDIF
	2052	write(,)trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
[5084]	2053	if (ncidpl.eq.-99) ncidpl=ncidQ
[1632]	2054	endif
	2055	! Pression de surface
[5084]	2056	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[3995]	2057	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5084]	2058	IF (rcode.NE.NF_NOERR) THEN
[3995]	2059	abort_message='Nudging: error -> no file ps.nc'
	2060	CALL abort_gcm(modname,abort_message,1)
	2061	ENDIF
	2062	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[5084]	2063	IF (rcode.NE.NF_NOERR) THEN
[3995]	2064	abort_message='Nudging: error -> no SP variable in file ps.nc'
	2065	CALL abort_gcm(modname,abort_message,1)
	2066	ENDIF
	2067	write(,)trim(modname)//' ncidps,varidps',ncidps,varidps
[1632]	2068	endif
	2069	! Coordonnee verticale
[5084]	2070	if (guide_plevs.EQ.0) then
[3995]	2071	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]	2072	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[3995]	2073	write(,)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1632]	2074	endif
	2075	! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]	2076	if (guide_plevs.EQ.1) then
	2077	#ifdef NC_DOUBLE
	2078	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
	2079	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
	2080	#else
	2081	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
	2082	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
	2083	#endif
	2084	elseif (guide_plevs.EQ.0) THEN
	2085	#ifdef NC_DOUBLE
	2086	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
	2087	#else
	2088	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
	2089	#endif
[1632]	2090	apnc=apnc*100.! conversion en Pascals
	2091	bpnc(:)=0.
	2092	endif
	2093	first=.FALSE.
	2094	endif ! (first)
	2095
	2096	! -----------------------------------------------------------------
	2097	! lecture des champs u, v, T, Q, ps
	2098	! -----------------------------------------------------------------
	2099
	2100	! dimensions pour les champs scalaires et le vent zonal
	2101	start(1)=1
	2102	start(2)=jjb_u
	2103	start(3)=1
	2104	start(4)=timestep
	2105
	2106	count(1)=1
	2107	count(2)=jjnb_u
	2108	count(3)=nlevnc
	2109	count(4)=1
	2110
	2111	IF (invert_y) start(2)=jjp1-jje_u+1
	2112	! Pression
[5084]	2113	if (guide_plevs.EQ.2) then
	2114	#ifdef NC_DOUBLE
	2115	status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,zu)
	2116	#else
	2117	status=NF_GET_VARA_REAL(ncidp,varidp,start,count,zu)
	2118	#endif
[1632]	2119	DO i=1,iip1
	2120	pnat2(i,:,:)=zu(:,:)
	2121	ENDDO
	2122
	2123	IF (invert_y) THEN
	2124	! PRINT*,"Invertion impossible actuellement"
	2125	! CALL abort_gcm(modname,abort_message,1)
	2126	CALL invert_lat(iip1,jjnb_u,nlevnc,pnat2)
	2127	ENDIF
	2128	endif
	2129	! Vent zonal
	2130	if (guide_u) then
[5084]	2131	#ifdef NC_DOUBLE
	2132	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
	2133	#else
	2134	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
	2135	#endif
[1632]	2136	DO i=1,iip1
	2137	unat2(i,:,:)=zu(:,:)
	2138	ENDDO
	2139
	2140	IF (invert_y) THEN
	2141	! PRINT*,"Invertion impossible actuellement"
	2142	! CALL abort_gcm(modname,abort_message,1)
	2143	CALL invert_lat(iip1,jjnb_u,nlevnc,unat2)
	2144	ENDIF
	2145	endif
	2146
[2028]	2147
[1632]	2148	! Temperature
	2149	if (guide_T) then
[5084]	2150	#ifdef NC_DOUBLE
	2151	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
	2152	#else
	2153	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
	2154	#endif
[1632]	2155	DO i=1,iip1
	2156	tnat2(i,:,:)=zu(:,:)
	2157	ENDDO
	2158
	2159	IF (invert_y) THEN
	2160	! PRINT*,"Invertion impossible actuellement"
	2161	! CALL abort_gcm(modname,abort_message,1)
	2162	CALL invert_lat(iip1,jjnb_u,nlevnc,tnat2)
	2163	ENDIF
	2164	endif
	2165
	2166	! Humidite
	2167	if (guide_Q) then
[5084]	2168	#ifdef NC_DOUBLE
	2169	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
	2170	#else
	2171	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
	2172	#endif
[1632]	2173	DO i=1,iip1
	2174	qnat2(i,:,:)=zu(:,:)
	2175	ENDDO
	2176
	2177	IF (invert_y) THEN
	2178	! PRINT*,"Invertion impossible actuellement"
	2179	! CALL abort_gcm(modname,abort_message,1)
	2180	CALL invert_lat(iip1,jjnb_u,nlevnc,qnat2)
	2181	ENDIF
	2182	endif
	2183
	2184	! Vent meridien
	2185	if (guide_v) then
	2186	start(2)=jjb_v
	2187	count(2)=jjnb_v
	2188	IF (invert_y) start(2)=jjm-jje_v+1
[5084]	2189	#ifdef NC_DOUBLE
	2190	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
	2191	#else
	2192	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
	2193	#endif
[1632]	2194	DO i=1,iip1
	2195	vnat2(i,:,:)=zv(:,:)
	2196	ENDDO
	2197
	2198	IF (invert_y) THEN
[2036]	2199
[1632]	2200	! PRINT*,"Invertion impossible actuellement"
	2201	! CALL abort_gcm(modname,abort_message,1)
	2202	CALL invert_lat(iip1,jjnb_v,nlevnc,vnat2)
	2203	ENDIF
	2204	endif
	2205
	2206	! Pression de surface
[5084]	2207	if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1632]	2208	start(2)=jjb_u
	2209	start(3)=timestep
	2210	start(4)=0
	2211	count(2)=jjnb_u
	2212	count(3)=1
	2213	count(4)=0
	2214	IF (invert_y) start(2)=jjp1-jje_u+1
[5084]	2215	#ifdef NC_DOUBLE
	2216	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
	2217	#else
	2218	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
	2219	#endif
[1632]	2220	DO i=1,iip1
	2221	psnat2(i,:)=zu(:,1)
	2222	ENDDO
	2223
	2224	IF (invert_y) THEN
	2225	! PRINT*,"Invertion impossible actuellement"
	2226	! CALL abort_gcm(modname,abort_message,1)
	2227	CALL invert_lat(iip1,jjnb_u,1,psnat2)
	2228	ENDIF
	2229	endif
	2230
	2231	END SUBROUTINE guide_read2D
	2232
	2233	!=======================================================================
[2028]	2234	SUBROUTINE guide_out(varname,hsize,vsize,field_loc,factt)
[1823]	2235	USE parallel_lmdz
[2028]	2236	USE mod_hallo, ONLY : gather_field_u, gather_field_v
[2597]	2237	USE comconst_mod, ONLY: pi
[2600]	2238	USE comvert_mod, ONLY: presnivs
[2740]	2239	use netcdf95, only: nf95_def_var, nf95_put_var
[5084]	2240	use netcdf, only: nf90_float
[2740]	2241
[1632]	2242	IMPLICIT NONE
	2243
	2244	INCLUDE "dimensions.h"
	2245	INCLUDE "paramet.h"
[5084]	2246	INCLUDE "netcdf.inc"
[1632]	2247	INCLUDE "comgeom2.h"
	2248
	2249	! Variables entree
[2028]	2250	CHARACTER(), INTENT(IN) :: varname
[1632]	2251	INTEGER, INTENT (IN) :: hsize,vsize
[2036]	2252	! REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field_loc
	2253	REAL, DIMENSION (:,:), INTENT(IN) :: field_loc
	2254	REAL factt
[1632]	2255
	2256	! Variables locales
	2257	INTEGER, SAVE :: timestep=0
	2258	! Identites fichier netcdf
	2259	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
	2260	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
[2740]	2261	INTEGER :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
[1632]	2262	INTEGER, DIMENSION (3) :: dim3
	2263	INTEGER, DIMENSION (4) :: dim4,count,start
[2028]	2264	INTEGER :: ierr, varid,l
[2740]	2265	REAL zu(ip1jmp1),zv(ip1jm), zt(iip1, jjp1), zq(iip1, jjp1)
[2036]	2266	REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: field_glo
[3995]	2267	CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
[1632]	2268
[2028]	2269	!$OMP MASTER
[2036]	2270	ALLOCATE(field_glo(iip1,hsize,vsize))
	2271	!$OMP END MASTER
	2272	!$OMP BARRIER
[2028]	2273
[3995]	2274	! write(,)trim(modname)//' after allocation ',hsize,vsize
[2036]	2275
[2028]	2276	IF (hsize==jjp1) THEN
[2036]	2277	CALL gather_field_u(field_loc,field_glo,vsize)
[2028]	2278	ELSE IF (hsize==jjm) THEN
	2279	CALL gather_field_v(field_loc,field_glo, vsize)
	2280	ENDIF
	2281
[3995]	2282	! write(,)trim(modname)//' after gather '
[2028]	2283	CALL Gather_field_u(alpha_u,zu,1)
[2740]	2284	CALL Gather_field_u(alpha_t,zt,1)
	2285	CALL Gather_field_u(alpha_q,zq,1)
[2028]	2286	CALL Gather_field_v(alpha_v,zv,1)
[2036]	2287
	2288	IF (mpi_rank > 0) THEN
	2289	!$OMP MASTER
	2290	DEALLOCATE(field_glo)
[2034]	2291	!$OMP END MASTER
	2292	!$OMP BARRIER
[2028]	2293
[2036]	2294	RETURN
	2295	ENDIF
[1632]	2296
[2034]	2297	!$OMP MASTER
[5084]	2298	IF (timestep.EQ.0) THEN
[1632]	2299	! ----------------------------------------------
	2300	! initialisation fichier de sortie
	2301	! ----------------------------------------------
	2302	! Ouverture du fichier
[3803]	2303	ierr=NF_CREATE("guide_ins.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
[1632]	2304	! Definition des dimensions
	2305	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
	2306	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
	2307	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
	2308	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
	2309	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
	2310	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
	2311
	2312	! Creation des variables dimensions
	2313	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
	2314	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
	2315	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
	2316	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
	2317	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
	2318	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
	2319	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
[2028]	2320	ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
	2321	ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
[2740]	2322	call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
	2323	varid_alpha_t)
	2324	call nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
	2325	varid_alpha_q)
	2326
[1632]	2327	ierr=NF_ENDDEF(nid)
	2328
	2329	! Enregistrement des variables dimensions
[5084]	2330	#ifdef NC_DOUBLE
	2331	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
	2332	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
	2333	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
	2334	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
	2335	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
	2336	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
	2337	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
	2338	ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,zu)
	2339	ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,zv)
	2340	#else
	2341	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
	2342	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
	2343	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
	2344	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
	2345	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
	2346	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
	2347	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
	2348	ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
	2349	ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
	2350	#endif
[2740]	2351	call nf95_put_var(nid, varid_alpha_t, zt)
	2352	call nf95_put_var(nid, varid_alpha_q, zq)
[1632]	2353	! --------------------------------------------------------------------
	2354	! Cr�ation des variables sauvegard�es
	2355	! --------------------------------------------------------------------
	2356	ierr = NF_REDEF(nid)
	2357	! Pressure (GCM)
	2358	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[2028]	2359	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
[1632]	2360	! Surface pressure (guidage)
	2361	IF (guide_P) THEN
	2362	dim3=(/id_lonv,id_latu,id_tim/)
	2363	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
	2364	ENDIF
	2365	! Zonal wind
	2366	IF (guide_u) THEN
	2367	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
[2028]	2368	ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
	2369	ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
[1632]	2370	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
	2371	ENDIF
	2372	! Merid. wind
	2373	IF (guide_v) THEN
	2374	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
[2028]	2375	ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
	2376	ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
[1632]	2377	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
	2378	ENDIF
	2379	! Pot. Temperature
	2380	IF (guide_T) THEN
	2381	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	2382	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
	2383	ENDIF
	2384	! Specific Humidity
	2385	IF (guide_Q) THEN
	2386	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
	2387	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
	2388	ENDIF
	2389
	2390	ierr = NF_ENDDEF(nid)
	2391	ierr = NF_CLOSE(nid)
	2392	ENDIF ! timestep=0
	2393
	2394	! --------------------------------------------------------------------
	2395	! Enregistrement du champ
	2396	! --------------------------------------------------------------------
	2397
	2398	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
	2399
[2028]	2400	IF (varname=="SP") timestep=timestep+1
	2401
	2402	ierr = NF_INQ_VARID(nid,varname,varid)
[1632]	2403	SELECT CASE (varname)
[2028]	2404	CASE ("SP","ps")
[1632]	2405	start=(/1,1,1,timestep/)
	2406	count=(/iip1,jjp1,llm,1/)
[2028]	2407	CASE ("v","va","vcov")
[1632]	2408	start=(/1,1,1,timestep/)
	2409	count=(/iip1,jjm,llm,1/)
[2028]	2410	CASE DEFAULT
[1632]	2411	start=(/1,1,1,timestep/)
	2412	count=(/iip1,jjp1,llm,1/)
[2028]	2413	END SELECT
	2414
[2036]	2415	!$OMP END MASTER
	2416	!$OMP BARRIER
	2417
[2028]	2418	SELECT CASE (varname)
	2419
	2420	CASE("u","ua")
[2036]	2421	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2039]	2422	DO l=1,llm
	2423	field_glo(:,2:jjm,l)=field_glo(:,2:jjm,l)/cu(:,2:jjm)
	2424	field_glo(:,1,l)=0. ; field_glo(:,jjp1,l)=0.
	2425	ENDDO
[2028]	2426	CASE("v","va")
[2036]	2427	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2039]	2428	DO l=1,llm
	2429	field_glo(:,:,l)=field_glo(:,:,l)/cv(:,:)
	2430	ENDDO
[2028]	2431	END SELECT
	2432
	2433	! if (varname=="ua") then
	2434	! call dump2d(iip1,jjp1,field_glo,'ua gui1 1ere couche ')
	2435	! call dump2d(iip1,jjp1,field_glo(:,:,llm),'ua gui1 llm ')
	2436	! endif
[2036]	2437
	2438	!$OMP MASTER
	2439
[5084]	2440	#ifdef NC_DOUBLE
	2441	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field_glo)
	2442	#else
	2443	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field_glo)
	2444	#endif
	2445
[1632]	2446	ierr = NF_CLOSE(nid)
	2447
[2036]	2448	DEALLOCATE(field_glo)
[2028]	2449	!$OMP END MASTER
	2450	!$OMP BARRIER
	2451
[1632]	2452	END SUBROUTINE guide_out
	2453
	2454
	2455	!===========================================================================
	2456	subroutine correctbid(iim,nl,x)
	2457	integer iim,nl
	2458	real x(iim+1,nl)
	2459	integer i,l
	2460	real zz
	2461
	2462	do l=1,nl
	2463	do i=2,iim-1
[5084]	2464	if(abs(x(i,l)).gt.1.e10) then
[1632]	2465	zz=0.5*(x(i-1,l)+x(i+1,l))
	2466	print*,'correction ',i,l,x(i,l),zz
	2467	x(i,l)=zz
	2468	endif
	2469	enddo
	2470	enddo
	2471	return
	2472	end subroutine correctbid
	2473
[2028]	2474
	2475	!====================================================================
	2476	! Ascii debug output. Could be reactivated
	2477	!====================================================================
	2478
	2479	subroutine dump2du(var,varname)
	2480	use parallel_lmdz
	2481	use mod_hallo
	2482	implicit none
	2483	include 'dimensions.h'
	2484	include 'paramet.h'
	2485
	2486	CHARACTER (len=*) :: varname
	2487
	2488
	2489	real, dimension(ijb_u:ije_u) :: var
	2490
	2491	real, dimension(ip1jmp1) :: var_glob
	2492
	2493	RETURN
	2494
	2495	call barrier
	2496	CALL Gather_field_u(var,var_glob,1)
	2497	call barrier
	2498
	2499	if (mpi_rank==0) then
	2500	call dump2d(iip1,jjp1,var_glob,varname)
	2501	endif
	2502
	2503	call barrier
	2504
	2505	return
	2506	end subroutine dump2du
	2507
	2508	!====================================================================
	2509	! Ascii debug output. Could be reactivated
	2510	!====================================================================
	2511	subroutine dumpall
	2512	implicit none
	2513	include "dimensions.h"
	2514	include "paramet.h"
	2515	include "comgeom.h"
	2516	call barrier
	2517	call dump2du(alpha_u(ijb_u:ije_u),' alpha_u couche 1')
	2518	call dump2du(unat2(:,jjbu:jjeu,nlevnc),' unat2 couche nlevnc')
	2519	call dump2du(ugui1(ijb_u:ije_u,1)*sqrt(unscu2(ijb_u:ije_u)),' ugui1 couche 1')
	2520	return
	2521	end subroutine dumpall
	2522
[1632]	2523	!===========================================================================
	2524	END MODULE guide_loc_mod

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

Download in other formats:

Original Format