Context Navigation

guide_mod.F90 @ 2041

Last change on this file since 2041 was 2025, checked in by fhourdin, 11 years ago
Enrichissement des sorties pour le guidage Enriched outputs for nudging
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 Property svn:executable set to ``*
File size: 58.0 KB

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

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

Download in other formats:

Original Format