Context Navigation

guide_loc_mod.F90 @ 5169

Last change on this file since 5169 was 5084, checked in by Laurent Fairhead, 12 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