Context Navigation

source: LMDZ5/branches/testing/libf/phylmd/1Dconv.h @ 2157

Last change on this file since 2157 was 2056, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1997:2055 into testing branch
Property svn:executable set to ``*
File size: 31.8 KB

Rev	Line
[2019]	1	subroutine get_uvd(itap,dtime,file_forctl,file_fordat, &
	2	& ht,hq,hw,hu,hv,hthturb,hqturb, &
	3	& Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
	4	!
[2017]	5	implicit none
	6
[2019]	7	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	8	! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
	9	! pouvoir calculer la convergence et le cisaillement dans la physiq
	10	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[2017]	11
	12	#include "YOMCST.h"
	13
	14	INTEGER klev
	15	REAL play(100) !pression en Pa au milieu de chaque couche GCM
	16	INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
	17	REAL coef1(100) !coefficient d interpolation
	18	REAL coef2(100) !coefficient d interpolation
	19
	20	INTEGER nblvlm !nombre de niveau de pression du mesoNH
	21	REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH
	22	REAL hplaym(100) !pression en hPa milieux des couches Meso-NH
	23
	24	integer i,j,k,ll,in
	25
	26	CHARACTER*80 file_forctl,file_fordat
	27
	28	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	29	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
	30
[2019]	31	!======================================================================
	32	! methode: on va chercher les donnees du mesoNH de meteo france, on y
	33	! a acces a tout pas detemps grace a la routine rdgrads qui
	34	! est une boucle lisant dans ces fichiers.
	35	! Puis on interpole ces donnes sur les 11 niveaux du gcm et
	36	! et sur les pas de temps de ce meme gcm
	37	!----------------------------------------------------------------------
	38	! input:
	39	! pasmax :nombre de pas de temps maximum du mesoNH
	40	! dt :pas de temps du meso_NH (en secondes)
	41	!----------------------------------------------------------------------
[2017]	42	integer pasmax,dt
	43	save pasmax,dt
[2019]	44	!----------------------------------------------------------------------
	45	! arguments:
	46	! itap :compteur de la physique(le nombre de ces pas est
	47	! fixe dans la subroutine calcul_ini_gcm de interpo
	48	! -lation
	49	! dtime :pas detemps du gcm (en secondes)
	50	! ht :convergence horizontale de temperature(K/s)
	51	! hq : " " d humidite (kg/kg/s)
	52	! hw :vitesse verticale moyenne (m/s**2)
	53	! hu :convergence horizontale d impulsion le long de x
	54	! (kg/(m^2 s^2)
	55	! hv : idem le long de y.
	56	! Ts : Temperature de surface (K)
	57	! imp_fcg: var. logical .eq. T si forcage en impulsion
	58	! ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
	59	! Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
	60	! Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
	61	!----------------------------------------------------------------------
[2017]	62	integer itap
	63	real dtime
	64	real ht(100)
	65	real hq(100)
	66	real hu(100)
	67	real hv(100)
	68	real hw(100)
	69	real hthturb(100)
	70	real hqturb(100)
	71	real Ts, Ts_subr
	72	logical imp_fcg
	73	logical ts_fcg
	74	logical Tp_fcg
	75	logical Turb_fcg
[2019]	76	!----------------------------------------------------------------------
	77	! Variables internes de get_uvd (note : l interpolation temporelle
	78	! est faite entre les pas de temps before et after, sur les variables
	79	! definies sur la grille du SCM; on atteint exactement les valeurs Meso
	80	! aux milieux des pas de temps Meso)
	81	! time0 :date initiale en secondes
	82	! time :temps associe a chaque pas du SCM
	83	! pas :numero du pas du meso_NH (on lit en pas : le premier pas
	84	! des donnees est duplique)
	85	! pasprev :numero du pas de lecture precedent
	86	! htaft :advection horizontale de temp. au pas de temps after
	87	! hqaft : " " d humidite "
	88	! hwaft :vitesse verticalle moyenne au pas de temps after
	89	! huaft,hvaft :advection horizontale d impulsion au pas de temps after
	90	! tsaft : surface temperature 'after time step'
	91	! htbef :idem htaft, mais pour le pas de temps before
	92	! hqbef :voir hqaft
	93	! hwbef :voir hwaft
	94	! hubef,hvbef : idem huaft,hvaft, mais pour before
	95	! tsbef : surface temperature 'before time step'
	96	!----------------------------------------------------------------------
[2017]	97	integer time0,pas,pasprev
	98	save time0,pas,pasprev
	99	real time
	100	real htaft(100),hqaft(100),hwaft(100),huaft(100),hvaft(100)
	101	real hthturbaft(100),hqturbaft(100)
	102	real Tsaft
	103	save htaft,hqaft,hwaft,huaft,hvaft,hthturbaft,hqturbaft
	104	real htbef(100),hqbef(100),hwbef(100),hubef(100),hvbef(100)
	105	real hthturbbef(100),hqturbbef(100)
	106	real Tsbef
	107	save htbef,hqbef,hwbef,hubef,hvbef,hthturbbef,hqturbbef
[2019]	108	!
[2017]	109	real timeaft,timebef
	110	save timeaft,timebef
	111	integer temps
	112	character*4 string
[2019]	113	!----------------------------------------------------------------------
	114	! variables arguments de la subroutine rdgrads
	115	!---------------------------------------------------------------------
[2017]	116	integer icompt,icomp1 !compteurs de rdgrads
	117	real z(100) ! altitude (grille Meso)
	118	real ht_mes(100) !convergence horizontale de temperature
	119	!-(grille Meso)
	120	real hq_mes(100) !convergence horizontale d humidite
	121	!(grille Meso)
	122	real hw_mes(100) !vitesse verticale moyenne
	123	!(grille Meso)
	124	real hu_mes(100),hv_mes(100) !convergence horizontale d impulsion
	125	!(grille Meso)
	126	real hthturb_mes(100) !tendance horizontale de T_pot, due aux
	127	!flux turbulents
	128	real hqturb_mes(100) !tendance horizontale d humidite, due aux
	129	!flux turbulents
[2019]	130	!
	131	!---------------------------------------------------------------------
	132	! variable argument de la subroutine copie
	133	!---------------------------------------------------------------------
	134	! SB real pplay(100) !pression en milieu de couche du gcm
	135	! SB !argument de la physique
	136	!---------------------------------------------------------------------
	137	! variables destinees a la lecture du pas de temps du fichier de donnees
	138	!---------------------------------------------------------------------
[2017]	139	character*80 aaa,atemps,spaces,apasmax
	140	integer nch,imn,ipa
[2019]	141	!---------------------------------------------------------------------
	142	! procedures appelees
[2017]	143	external rdgrads !lire en iterant dans forcing.dat
[2019]	144	!---------------------------------------------------------------------
[2017]	145	print*,'le pas itap est:',itap
[2019]	146	!* on determine le pas du meso_NH correspondant au nouvel itap *
	147	!* pour aller chercher les champs dans rdgrads *
	148	!
[2017]	149	time=time0+itap*dtime
[2019]	150	!c temps=int(time/dt+1)
	151	!c pas=min(temps,pasmax)
[2017]	152	temps = 1 + int((dt + 2time)/(2dt))
	153	pas=min(temps,pasmax-1)
	154	print*,'le pas Meso est:',pas
[2019]	155	!
	156	!
	157	!===================================================================
	158	!
	159	!* on remplit les champs before avec les champs after du pas *
	160	!* precedent en format gcm *
[2017]	161	if(pas.gt.pasprev)then
	162	do i=1,klev
	163	htbef(i)=htaft(i)
	164	hqbef(i)=hqaft(i)
	165	hwbef(i)=hwaft(i)
	166	hubef(i)=huaft(i)
	167	hvbef(i)=hvaft(i)
	168	hThTurbbef(i)=hThTurbaft(i)
	169	hqTurbbef(i)=hqTurbaft(i)
	170	enddo
	171	tsbef = tsaft
	172	timebef=pasprev*dt
	173	timeaft=timebef+dt
	174	icomp1 = nblvlm*4
	175	IF (ts_fcg) icomp1 = icomp1 + 1
	176	IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
	177	IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
	178	icompt = icomp1*pas
	179	print *, 'imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt'
	180	print *, imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt
	181	print*,'le pas pas est:',pas
[2019]	182	!* on va chercher les nouveaux champs after dans toga.dat *
	183	!* champs en format meso_NH *
	184	open(99,FILE=file_fordat,FORM='UNFORMATTED', &
	185	& ACCESS='DIRECT',RECL=8)
	186	call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes &
	187	& ,hu_mes,hv_mes,hthturb_mes,hqturb_mes &
	188	& ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
	189	!
[2017]	190
	191	if(Tp_fcg) then
[2019]	192	! (le forcage est donne en temperature potentielle)
[2017]	193	do i = 1,nblvlm
	194	ht_mes(i) = ht_mes(i)(hplaym(i)/1000.)*rkappa
	195	enddo
	196	endif ! Tp_fcg
	197	if(Turb_fcg) then
	198	do i = 1,nblvlm
	199	hThTurb_mes(i) = hThTurb_mes(i)(hplaym(i)/1000.)*rkappa
	200	enddo
	201	endif ! Turb_fcg
[2019]	202	!
[2017]	203	print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
	204	print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
	205	print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
	206	if(imp_fcg) then
	207	print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
	208	print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
	209	endif
	210	if(Turb_fcg) then
	211	print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
	212	print*,'hqTurb_mes ',(hqTurb_mes(i),i=1,nblvlm)
	213	endif
	214	IF (ts_fcg) print*,'ts_subr', ts_subr
[2019]	215	!* on interpole les champs meso_NH sur les niveaux de pression*
	216	!* gcm . on obtient le nouveau champ after *
[2017]	217	do k=1,klev
	218	if (JM(k) .eq. 0) then
	219	htaft(k)= ht_mes(jm(k)+1)
	220	hqaft(k)= hq_mes(jm(k)+1)
	221	hwaft(k)= hw_mes(jm(k)+1)
	222	if(imp_fcg) then
	223	huaft(k)= hu_mes(jm(k)+1)
	224	hvaft(k)= hv_mes(jm(k)+1)
	225	endif ! imp_fcg
	226	if(Turb_fcg) then
	227	hThTurbaft(k)= hThTurb_mes(jm(k)+1)
	228	hqTurbaft(k)= hqTurb_mes(jm(k)+1)
	229	endif ! Turb_fcg
	230	else ! JM(k) .eq. 0
	231	htaft(k)=coef1(k)ht_mes(jm(k))+coef2(k)ht_mes(jm(k)+1)
	232	hqaft(k)=coef1(k)hq_mes(jm(k))+coef2(k)hq_mes(jm(k)+1)
	233	hwaft(k)=coef1(k)hw_mes(jm(k))+coef2(k)hw_mes(jm(k)+1)
	234	if(imp_fcg) then
	235	huaft(k)=coef1(k)hu_mes(jm(k))+coef2(k)hu_mes(jm(k)+1)
	236	hvaft(k)=coef1(k)hv_mes(jm(k))+coef2(k)hv_mes(jm(k)+1)
	237	endif ! imp_fcg
	238	if(Turb_fcg) then
[2019]	239	hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k)) &
	240	& +coef2(k)*hThTurb_mes(jm(k)+1)
	241	hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k)) &
	242	& +coef2(k)*hqTurb_mes(jm(k)+1)
[2017]	243	endif ! Turb_fcg
	244	endif ! JM(k) .eq. 0
	245	enddo
	246	tsaft = ts_subr
	247	pasprev=pas
	248	else ! pas.gt.pasprev
	249	print*,'timebef est:',timebef
	250	endif ! pas.gt.pasprev fin du bloc relatif au passage au pas
	251	!de temps (meso) suivant
[2019]	252	!* si on atteint le pas max des donnees experimentales ,on *
	253	!* on conserve les derniers champs calcules *
[2017]	254	if(temps.ge.pasmax)then
	255	do ll=1,klev
	256	ht(ll)=htaft(ll)
	257	hq(ll)=hqaft(ll)
	258	hw(ll)=hwaft(ll)
	259	hu(ll)=huaft(ll)
	260	hv(ll)=hvaft(ll)
	261	hThTurb(ll)=hThTurbaft(ll)
	262	hqTurb(ll)=hqTurbaft(ll)
	263	enddo
	264	ts_subr = tsaft
	265	else ! temps.ge.pasmax
[2019]	266	!* on interpole sur les pas de temps de 10mn du gcm a partir *
	267	! des pas de temps de 1h du meso_NH *
[2017]	268	do j=1,klev
	269	ht(j)=((timeaft-time)htbef(j)+(time-timebef)htaft(j))/dt
	270	hq(j)=((timeaft-time)hqbef(j)+(time-timebef)hqaft(j))/dt
	271	hw(j)=((timeaft-time)hwbef(j)+(time-timebef)hwaft(j))/dt
	272	if(imp_fcg) then
	273	hu(j)=((timeaft-time)hubef(j)+(time-timebef)huaft(j))/dt
	274	hv(j)=((timeaft-time)hvbef(j)+(time-timebef)hvaft(j))/dt
	275	endif ! imp_fcg
	276	if(Turb_fcg) then
[2019]	277	hThTurb(j)=((timeaft-time)*hThTurbbef(j) &
	278	& +(time-timebef)*hThTurbaft(j))/dt
	279	hqTurb(j)= ((timeaft-time)*hqTurbbef(j) &
	280	& +(time-timebef)*hqTurbaft(j))/dt
[2017]	281	endif ! Turb_fcg
	282	enddo
	283	ts_subr = ((timeaft-time)tsbef + (time-timebef)tsaft)/dt
	284	endif ! temps.ge.pasmax
[2019]	285	!
[2017]	286	print *,' time,timebef,timeaft',time,timebef,timeaft
	287	print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
	288	do j= 1,klev
[2019]	289	print *, j,ht(j),htbef(j),htaft(j), &
	290	& hthturb(j),hthturbbef(j),hthturbaft(j)
[2017]	291	enddo
	292	print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
	293	do j= 1,klev
[2019]	294	print *, j,hq(j),hqbef(j),hqaft(j), &
	295	& hqturb(j),hqturbbef(j),hqturbaft(j)
[2017]	296	enddo
[2019]	297	!
	298	!-------------------------------------------------------------------
	299	!
[2017]	300	IF (Ts_fcg) Ts = Ts_subr
	301	return
[2019]	302	!
	303	!-----------------------------------------------------------------------
	304	! on sort les champs de "convergence" pour l instant initial 'in'
	305	! ceci se passe au pas temps itap=0 de la physique
	306	!-----------------------------------------------------------------------
	307	entry get_uvd2(itap,dtime,file_forctl,file_fordat, &
	308	& ht,hq,hw,hu,hv,hThTurb,hqTurb,ts, &
	309	& imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
[2017]	310	print*,'le pas itap est:',itap
[2019]	311	!
	312	!===================================================================
	313	!
[2017]	314	write(*,'(a)') 'OPEN '//file_forctl
	315	open(97,FILE=file_forctl,FORM='FORMATTED')
[2019]	316	!
	317	!------------------
[2017]	318	do i=1,1000
	319	read(97,1000,end=999) string
	320	1000 format (a4)
	321	if (string .eq. 'TDEF') go to 50
	322	enddo
	323	50 backspace(97)
[2019]	324	!-------------------------------------------------------------------
	325	! * on lit le pas de temps dans le fichier de donnees *
	326	! * "forcing.ctl" et pasmax *
	327	!-------------------------------------------------------------------
[2017]	328	read(97,2000) aaa
	329	2000 format (a80)
	330	print*,'aaa est',aaa
	331	aaa=spaces(aaa,1)
	332	print*,'aaa',aaa
	333	call getsch(aaa,' ',' ',5,atemps,nch)
	334	print*,'atemps est',atemps
	335	atemps=atemps(1:nch-2)
	336	print*,'atemps',atemps
	337	read(atemps,*) imn
	338	dt=imn*60
	339	print*,'le pas de temps dt',dt
	340	call getsch(aaa,' ',' ',2,apasmax,nch)
	341	apasmax=apasmax(1:nch)
	342	read(apasmax,*) ipa
	343	pasmax=ipa
	344	print*,'pasmax est',pasmax
	345	CLOSE(97)
[2019]	346	!------------------------------------------------------------------
	347	! * on lit le pas de temps initial de la simulation *
	348	!------------------------------------------------------------------
[2017]	349	in=itap
[2019]	350	!c pasprev=in
	351	!c time0=dt*(pasprev-1)
[2017]	352	pasprev=in-1
	353	time0=dt*pasprev
[2019]	354	!
[2017]	355	close(98)
[2019]	356	!
[2017]	357	write(*,'(a)') 'OPEN '//file_fordat
[2019]	358	open(99,FILE=file_fordat,FORM='UNFORMATTED', &
	359	& ACCESS='DIRECT',RECL=8)
[2017]	360	icomp1 = nblvlm*4
	361	IF (ts_fcg) icomp1 = icomp1 + 1
	362	IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
	363	IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
	364	icompt = icomp1*(in-1)
[2019]	365	call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes &
	366	& ,hu_mes,hv_mes,hthturb_mes,hqturb_mes &
	367	& ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
[2017]	368	print *, 'get_uvd : rdgrads ->'
	369	print *, tp_fcg
[2019]	370	!
	371	! following commented out because we have temperature already in ARM case
	372	! (otherwise this is the potential temperature )
	373	!------------------------------------------------------------------------
[2017]	374	if(Tp_fcg) then
	375	do i = 1,nblvlm
	376	ht_mes(i) = ht_mes(i)(hplaym(i)/1000.)*rkappa
	377	enddo
	378	endif ! Tp_fcg
	379	print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
	380	print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
	381	print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
	382	if(imp_fcg) then
	383	print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
	384	print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
	385	print*,'t',ts_subr
	386	endif ! imp_fcg
	387	if(Turb_fcg) then
	388	print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
	389	print*,'hqTurb ', (hqTurb_mes(i),i=1,nblvlm)
	390	endif ! Turb_fcg
[2019]	391	!----------------------------------------------------------------------
	392	! on a obtenu des champs initiaux sur les niveaux du meso_NH
	393	! on interpole sur les niveaux du gcm(niveau pression bien sur!)
	394	!-----------------------------------------------------------------------
[2017]	395	do k=1,klev
	396	if (JM(k) .eq. 0) then
[2019]	397	!FKC bug? ne faut il pas convertir tsol en tendance ????
	398	!RT bug taken care of by removing the stuff
[2017]	399	htaft(k)= ht_mes(jm(k)+1)
	400	hqaft(k)= hq_mes(jm(k)+1)
	401	hwaft(k)= hw_mes(jm(k)+1)
	402	if(imp_fcg) then
	403	huaft(k)= hu_mes(jm(k)+1)
	404	hvaft(k)= hv_mes(jm(k)+1)
	405	endif ! imp_fcg
	406	if(Turb_fcg) then
	407	hThTurbaft(k)= hThTurb_mes(jm(k)+1)
	408	hqTurbaft(k)= hqTurb_mes(jm(k)+1)
	409	endif ! Turb_fcg
	410	else ! JM(k) .eq. 0
	411	htaft(k)=coef1(k)ht_mes(jm(k))+coef2(k)ht_mes(jm(k)+1)
	412	hqaft(k)=coef1(k)hq_mes(jm(k))+coef2(k)hq_mes(jm(k)+1)
	413	hwaft(k)=coef1(k)hw_mes(jm(k))+coef2(k)hw_mes(jm(k)+1)
	414	if(imp_fcg) then
	415	huaft(k)=coef1(k)hu_mes(jm(k))+coef2(k)hu_mes(jm(k)+1)
	416	hvaft(k)=coef1(k)hv_mes(jm(k))+coef2(k)hv_mes(jm(k)+1)
	417	endif ! imp_fcg
	418	if(Turb_fcg) then
[2019]	419	hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k)) &
	420	& +coef2(k)*hThTurb_mes(jm(k)+1)
	421	hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k)) &
	422	& +coef2(k)*hqTurb_mes(jm(k)+1)
[2017]	423	endif ! Turb_fcg
	424	endif ! JM(k) .eq. 0
	425	enddo
	426	tsaft = ts_subr
[2019]	427	! valeurs initiales des champs de convergence
[2017]	428	do k=1,klev
	429	ht(k)=htaft(k)
	430	hq(k)=hqaft(k)
	431	hw(k)=hwaft(k)
	432	if(imp_fcg) then
	433	hu(k)=huaft(k)
	434	hv(k)=hvaft(k)
	435	endif ! imp_fcg
	436	if(Turb_fcg) then
	437	hThTurb(k)=hThTurbaft(k)
	438	hqTurb(k) =hqTurbaft(k)
	439	endif ! Turb_fcg
	440	enddo
	441	ts_subr = tsaft
	442	close(99)
	443	close(98)
[2019]	444	!
	445	!-------------------------------------------------------------------
	446	!
	447	!
[2017]	448	100 IF (Ts_fcg) Ts = Ts_subr
	449	return
[2019]	450	!
[2017]	451	999 continue
	452	stop 'erreur lecture, file forcing.ctl'
	453	end
	454
[2019]	455	SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f &
	456	& ,d_t_adv,d_q_adv)
[2017]	457	use dimphy
	458	implicit none
	459
	460	#include "dimensions.h"
[2019]	461	!cccc#include "dimphy.h"
[2017]	462
	463	integer k
	464	real dtime, fact, du, dv, cx, cy, alx, aly
	465	real zt(klev), zq(klev,3)
[2019]	466	real vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
[2017]	467
	468	real d_t_adv(klev), d_q_adv(klev,3)
	469
[2019]	470	! Velocity of moving cell
[2017]	471	data cx,cy /12., -2./
	472
[2019]	473	! Dimensions of moving cell
	474	data alx,aly /100000.,150000./
[2017]	475
	476	do k = 1, klev
	477	du = abs(vu_f(k)-cx)/alx
	478	dv = abs(vv_f(k)-cy)/aly
	479	fact = dtime (du+dv-dudv*dtime)
	480	d_t_adv(k) = fact * (t_f(k)-zt(k))
	481	d_q_adv(k,1) = fact * (q_f(k,1)-zq(k,1))
	482	d_q_adv(k,2) = fact * (q_f(k,2)-zq(k,2))
	483	d_q_adv(k,3) = fact * (q_f(k,3)-zq(k,3))
	484	enddo
	485
	486	return
	487	end
	488
	489	SUBROUTINE copie(klevgcm,playgcm,psolgcm,file_forctl)
	490	implicit none
	491
[2019]	492	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	493	! cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
	494	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[2017]	495
	496	INTEGER klev !nombre de niveau de pression du GCM
	497	REAL play(100) !pression en Pa au milieu de chaque couche GCM
	498	INTEGER JM(100)
	499	REAL coef1(100) !coefficient d interpolation
	500	REAL coef2(100) !coefficient d interpolation
	501
	502	INTEGER nblvlm !nombre de niveau de pression du mesoNH
	503	REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH
	504	REAL hplaym(100)!pression en hecto-Pa des milieux de couche Meso-NH
	505
	506	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	507	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
	508
	509	integer k,klevgcm
	510	real playgcm(klevgcm) ! pression en milieu de couche du gcm
	511	real psolgcm
	512	character*80 file_forctl
	513
	514	klev = klevgcm
	515
[2019]	516	!---------------------------------------------------------------------
	517	! pression au milieu des couches du gcm dans la physiq
	518	! (SB: remplace le call conv_lipress_gcm(playgcm) )
	519	!---------------------------------------------------------------------
[2017]	520
	521	do k = 1, klev
	522	play(k) = playgcm(k)
	523	print*,'la pression gcm est:',play(k)
	524	enddo
	525
[2019]	526	!----------------------------------------------------------------------
	527	! lecture du descripteur des donnees Meso-NH (forcing.ctl):
	528	! -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
	529	! (on remplit le COMMON com2_phys_gcss)
	530	!----------------------------------------------------------------------
[2017]	531
	532	call mesolupbis(file_forctl)
	533
	534	print*,'la valeur de nblvlm est:',nblvlm
	535
[2019]	536	!----------------------------------------------------------------------
	537	! etude de la correspondance entre les niveaux meso.NH et GCM;
	538	! calcul des coefficients d interpolation coef1 et coef2
	539	! (on remplit le COMMON com1_phys_gcss)
	540	!----------------------------------------------------------------------
[2017]	541
	542	call corresbis(psolgcm)
	543
[2019]	544	!---------------------------------------------------------
	545	! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
	546	!---------------------------------------------------------
[2017]	547
	548	write(,) ' '
	549	write(,) 'TESTS com1_phys_gcss et com2_phys_gcss dans copie.F'
	550	write(,) '--------------------------------------'
	551	write(,) 'GCM: nb niveaux:',klev,' et pression, coeffs:'
	552	do k = 1, klev
	553	write(,) play(k), coef1(k), coef2(k)
	554	enddo
	555	write(,) 'MESO-NH: nb niveaux:',nblvlm,' et pression:'
	556	do k = 1, nblvlm
	557	write(,) playm(k), hplaym(k)
	558	enddo
	559	write(,) ' '
	560
	561	end
	562	SUBROUTINE mesolupbis(file_forctl)
	563	implicit none
[2019]	564	!
	565	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	566	!
	567	! Lecture descripteur des donnees MESO-NH (forcing.ctl):
	568	! -------------------------------------------------------
	569	!
	570	! Cette subroutine lit dans le fichier de controle "essai.ctl"
	571	! et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
	572	! des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
	573	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	574	!
[2017]	575	INTEGER nblvlm !nombre de niveau de pression du mesoNH
	576	REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH
	577	REAL hplaym(100) !pression en hPa des milieux de couche Meso-NH
	578	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
	579
	580	INTEGER i,lu,mlz,mlzh
	581
	582	character*80 file_forctl
	583
	584	character*4 a
	585	character*80 aaa,anblvl,spaces
	586	integer nch
	587
	588	lu=9
	589	open(lu,file=file_forctl,form='formatted')
[2019]	590	!
[2017]	591	do i=1,1000
	592	read(lu,1000,end=999) a
	593	if (a .eq. 'ZDEF') go to 100
	594	enddo
[2019]	595	!
[2017]	596	100 backspace(lu)
	597	print*,' DESCRIPTION DES 2 MODELES : '
	598	print*,' '
[2019]	599	!
[2017]	600	read(lu,2000) aaa
	601	2000 format (a80)
	602	aaa=spaces(aaa,1)
	603	call getsch(aaa,' ',' ',2,anblvl,nch)
	604	read(anblvl,*) nblvlm
	605
[2019]	606	!
[2017]	607	print*,'nbre de niveaux de pression Meso-NH :',nblvlm
	608	print*,' '
	609	print*,'pression en Pa de chaque couche du meso-NH :'
[2019]	610	!
[2017]	611	read(lu,*) (playm(mlz),mlz=1,nblvlm)
[2019]	612	! Si la pression est en HPa, la multiplier par 100
[2017]	613	if (playm(1) .lt. 10000.) then
	614	do mlz = 1,nblvlm
	615	playm(mlz) = playm(mlz)*100.
	616	enddo
	617	endif
	618	print*,(playm(mlz),mlz=1,nblvlm)
[2019]	619	!
[2017]	620	1000 format (a4)
	621	1001 format(5x,i2)
[2019]	622	!
[2017]	623	print*,' '
	624	do mlzh=1,nblvlm
	625	hplaym(mlzh)=playm(mlzh)/100.
	626	enddo
[2019]	627	!
[2017]	628	print*,'pression en hPa de chaque couche du meso-NH: '
	629	print*,(hplaym(mlzh),mlzh=1,nblvlm)
[2019]	630	!
[2017]	631	close (lu)
	632	return
[2019]	633	!
[2017]	634	999 stop 'erreur lecture des niveaux pression des donnees'
	635	end
	636
[2019]	637	SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur, &
	638	& ts_fcg,ts,imp_fcg,Turb_fcg)
[2017]	639	IMPLICIT none
	640	INTEGER itape,icount,icomp, nl
	641	real z(nl),ht(nl),hq(nl),hw(nl),hu(nl),hv(nl)
	642	real hthtur(nl),hqtur(nl)
	643	real ts
[2019]	644	!
[2017]	645	INTEGER k
[2019]	646	!
[2017]	647	LOGICAL imp_fcg,ts_fcg,Turb_fcg
[2019]	648	!
[2017]	649	icomp = icount
[2019]	650	!
	651	!
[2017]	652	do k=1,nl
	653	icomp=icomp+1
	654	read(itape,rec=icomp)z(k)
	655	print *,'icomp,k,z(k) ',icomp,k,z(k)
	656	enddo
	657	do k=1,nl
	658	icomp=icomp+1
	659	read(itape,rec=icomp)hT(k)
	660	print*, hT(k), k
	661	enddo
	662	do k=1,nl
	663	icomp=icomp+1
	664	read(itape,rec=icomp)hQ(k)
	665	enddo
[2019]	666	!
[2017]	667	if(turb_fcg) then
	668	do k=1,nl
	669	icomp=icomp+1
	670	read(itape,rec=icomp)hThTur(k)
	671	enddo
	672	do k=1,nl
	673	icomp=icomp+1
	674	read(itape,rec=icomp)hqTur(k)
	675	enddo
	676	endif
	677	print *,' apres lecture hthtur, hqtur'
[2019]	678	!
[2017]	679	if(imp_fcg) then
	680
	681	do k=1,nl
	682	icomp=icomp+1
	683	read(itape,rec=icomp)hu(k)
	684	enddo
	685	do k=1,nl
	686	icomp=icomp+1
	687	read(itape,rec=icomp)hv(k)
	688	enddo
	689
	690	endif
[2019]	691	!
[2017]	692	do k=1,nl
	693	icomp=icomp+1
	694	read(itape,rec=icomp)hw(k)
	695	enddo
[2019]	696	!
[2017]	697	if(ts_fcg) then
	698	icomp=icomp+1
	699	read(itape,rec=icomp)ts
	700	endif
[2019]	701	!
[2017]	702	print *,' rdgrads ->'
	703
	704	RETURN
	705	END
	706
	707	SUBROUTINE corresbis(psol)
	708	implicit none
	709
[2019]	710	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	711	! Cette subroutine calcule et affiche les valeurs des coefficients
	712	! d interpolation qui serviront dans la formule d interpolation elle-
	713	! meme.
	714	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[2017]	715
	716	INTEGER klev !nombre de niveau de pression du GCM
	717	REAL play(100) !pression en Pa au milieu de chaque couche GCM
	718	INTEGER JM(100)
	719	REAL coef1(100) !coefficient d interpolation
	720	REAL coef2(100) !coefficient d interpolation
	721
	722	INTEGER nblvlm !nombre de niveau de pression du mesoNH
	723	REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH
	724	REAL hplaym(100)!pression en hPa des milieux de couche Meso-NH
	725
	726	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	727	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
	728
	729	REAL psol
	730	REAL val
	731	INTEGER k, mlz
	732
	733
	734	do k=1,klev
	735	val=play(k)
	736	if (val .gt. playm(1)) then
	737	mlz = 0
	738	JM(1) = mlz
[2019]	739	coef1(1)=(playm(mlz+1)-val)/(playm(mlz+1)-psol)
	740	coef2(1)=(val-psol)/(playm(mlz+1)-psol)
[2017]	741	else if (val .gt. playm(nblvlm)) then
	742	do mlz=1,nblvlm
[2019]	743	if ( val .le. playm(mlz).and. val .gt. playm(mlz+1))then
[2017]	744	JM(k)=mlz
[2019]	745	coef1(k)=(playm(mlz+1)-val)/(playm(mlz+1)-playm(mlz))
	746	coef2(k)=(val-playm(mlz))/(playm(mlz+1)-playm(mlz))
[2017]	747	endif
	748	enddo
	749	else
	750	JM(k) = nblvlm-1
	751	coef1(k) = 0.
	752	coef2(k) = 0.
	753	endif
	754	enddo
[2019]	755	!
	756	!c if (play(klev) .le. playm(nblvlm)) then
	757	!c mlz=nblvlm-1
	758	!c JM(klev)=mlz
	759	!c coef1(klev)=(playm(mlz+1)-val)
	760	!c * /(playm(mlz+1)-playm(mlz))
	761	!c coef2(klev)=(val-playm(mlz))
	762	!c * /(playm(mlz+1)-playm(mlz))
	763	!c endif
	764	!
[2017]	765	print*,' '
	766	print*,' INTERPOLATION : '
	767	print*,' '
	768	print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
	769	print*,(JM(k),k=1,klev)
	770	print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
	771	print*,(JM(k),k=1,klev)
	772	print*,' '
[2019]	773	print*,'vals du premier coef d"interpolation pour les 9 niveaux: '
[2017]	774	print*,(coef1(k),k=1,klev)
	775	print*,' '
[2019]	776	print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:'
[2017]	777	print*,(coef2(k),k=1,klev)
[2019]	778	!
[2017]	779	return
	780	end
	781	SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
[2019]	782	!***************************************************************
	783	!* *
	784	!* *
	785	!* GETSCH *
	786	!* *
	787	!* *
	788	!* modified by : *
	789	!***************************************************************
	790	!* Return in SST the character string found between the NTH-1 and NTH
	791	!* occurence of the delimiter 'DEL' but before the terminator 'TRM' in
	792	!* the input string 'STR'. If TRM=DEL then STR is considered unlimited.
	793	!* NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
	794	!* NTH is greater than the number of delimiters in STR.
[2017]	795	IMPLICIT INTEGER (A-Z)
	796	CHARACTER STR(),DEL1,TRM1,SST()
	797	NCH=-1
	798	SST=' '
	799	IF(NTH.GT.0) THEN
	800	IF(TRM.EQ.DEL) THEN
	801	LENGTH=LEN(STR)
	802	ELSE
	803	LENGTH=INDEX(STR,TRM)-1
	804	IF(LENGTH.LT.0) LENGTH=LEN(STR)
	805	ENDIF
[2019]	806	!* Find beginning and end of the NTH DEL-limited substring in STR
[2017]	807	END=-1
	808	DO 1,N=1,NTH
	809	IF(END.EQ.LENGTH) RETURN
	810	BEG=END+2
	811	END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
	812	IF(END.EQ.BEG-2) END=LENGTH
[2019]	813	!* PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
[2017]	814	1 CONTINUE
	815	NCH=END-BEG+1
	816	IF(NCH.GT.0) SST=STR(BEG:END)
	817	ENDIF
	818	END
	819	CHARACTER() FUNCTION SPACES(STR,NSPACE)
[2019]	820	!
	821	! CERN PROGLIB# M433 SPACES .VERSION KERNFOR 4.14 860211
	822	! ORIG. 6/05/86 M.GOOSSENS/DD
	823	!
	824	!- The function value SPACES returns the character string STR with
	825	!- leading blanks removed and each occurence of one or more blanks
	826	!- replaced by NSPACE blanks inside the string STR
	827	!
[2017]	828	CHARACTER() STR
[2019]	829	!
[2017]	830	LENSPA = LEN(SPACES)
	831	SPACES = ' '
	832	IF (NSPACE.LT.0) NSPACE = 0
	833	IBLANK = 1
	834	ISPACE = 1
	835	100 INONBL = INDEXC(STR(IBLANK:),' ')
	836	IF (INONBL.EQ.0) THEN
	837	SPACES(ISPACE:) = STR(IBLANK:)
	838	GO TO 999
	839	ENDIF
	840	INONBL = INONBL + IBLANK - 1
	841	IBLANK = INDEX(STR(INONBL:),' ')
	842	IF (IBLANK.EQ.0) THEN
	843	SPACES(ISPACE:) = STR(INONBL:)
	844	GO TO 999
	845	ENDIF
	846	IBLANK = IBLANK + INONBL - 1
	847	SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
	848	ISPACE = ISPACE + IBLANK - INONBL + NSPACE
	849	IF (ISPACE.LE.LENSPA) GO TO 100
	850	999 END
	851	FUNCTION INDEXC(STR,SSTR)
[2019]	852	!
	853	! CERN PROGLIB# M433 INDEXC .VERSION KERNFOR 4.14 860211
	854	! ORIG. 26/03/86 M.GOOSSENS/DD
	855	!
	856	!- Find the leftmost position where substring SSTR does not match
	857	!- string STR scanning forward
	858	!
[2017]	859	CHARACTER() STR,SSTR
[2019]	860	!
[2017]	861	LENS = LEN(STR)
	862	LENSS = LEN(SSTR)
[2019]	863	!
[2017]	864	DO 10 I=1,LENS-LENSS+1
	865	IF (STR(I:I+LENSS-1).NE.SSTR) THEN
	866	INDEXC = I
	867	GO TO 999
	868	ENDIF
	869	10 CONTINUE
	870	INDEXC = 0
[2019]	871	!
[2017]	872	999 END

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

Download in other formats: