Context Navigation

source: LMDZ5/trunk/libf/phy1d/1Dconv.h @ 1963

Last change on this file since 1963 was 1960, checked in by fhourdin, 11 years ago

Nettoyage du code 1D pour limiter les warning en mode debug.
Essentiellement des declarations inutiles.

Cleaning of the 1D code in order to limit warnings in debug mode.
Essentially supressing unused declared variables.

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: 31.4 KB

Rev	Line
[1607]	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	c
	5	implicit none
	6
	7	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	8	c cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
	9	c pouvoir calculer la convergence et le cisaillement dans la physiq
	10	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	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
[1960]	24	integer i,j,k,ll,in
[1607]	25
	26	CHARACTER*80 file_forctl,file_fordat
	27
[1960]	28	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	29	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
[1607]	30
	31	c======================================================================
	32	c methode: on va chercher les donnees du mesoNH de meteo france, on y
	33	c a acces a tout pas detemps grace a la routine rdgrads qui
	34	c est une boucle lisant dans ces fichiers.
	35	c Puis on interpole ces donnes sur les 11 niveaux du gcm et
	36	c et sur les pas de temps de ce meme gcm
	37	c----------------------------------------------------------------------
	38	c input:
	39	c pasmax :nombre de pas de temps maximum du mesoNH
	40	c dt :pas de temps du meso_NH (en secondes)
	41	c----------------------------------------------------------------------
	42	integer pasmax,dt
	43	save pasmax,dt
	44	c----------------------------------------------------------------------
	45	c arguments:
	46	c itap :compteur de la physique(le nombre de ces pas est
	47	c fixe dans la subroutine calcul_ini_gcm de interpo
	48	c -lation
	49	c dtime :pas detemps du gcm (en secondes)
	50	c ht :convergence horizontale de temperature(K/s)
	51	c hq : " " d humidite (kg/kg/s)
	52	c hw :vitesse verticale moyenne (m/s**2)
	53	c hu :convergence horizontale d impulsion le long de x
	54	c (kg/(m^2 s^2)
	55	c hv : idem le long de y.
	56	c Ts : Temperature de surface (K)
	57	c imp_fcg: var. logical .eq. T si forcage en impulsion
	58	c ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
	59	c Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
	60	c Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
	61	c----------------------------------------------------------------------
	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
	76	c----------------------------------------------------------------------
	77	c Variables internes de get_uvd (note : l interpolation temporelle
	78	c est faite entre les pas de temps before et after, sur les variables
	79	c definies sur la grille du SCM; on atteint exactement les valeurs Meso
	80	c aux milieux des pas de temps Meso)
	81	c time0 :date initiale en secondes
	82	c time :temps associe a chaque pas du SCM
	83	c pas :numero du pas du meso_NH (on lit en pas : le premier pas
	84	c des donnees est duplique)
	85	c pasprev :numero du pas de lecture precedent
	86	c htaft :advection horizontale de temp. au pas de temps after
	87	c hqaft : " " d humidite "
	88	c hwaft :vitesse verticalle moyenne au pas de temps after
	89	c huaft,hvaft :advection horizontale d impulsion au pas de temps after
	90	c tsaft : surface temperature 'after time step'
	91	c htbef :idem htaft, mais pour le pas de temps before
	92	c hqbef :voir hqaft
	93	c hwbef :voir hwaft
	94	c hubef,hvbef : idem huaft,hvaft, mais pour before
	95	c tsbef : surface temperature 'before time step'
	96	c----------------------------------------------------------------------
	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
	108	c
	109	real timeaft,timebef
	110	save timeaft,timebef
	111	integer temps
	112	character*4 string
	113	c----------------------------------------------------------------------
	114	c variables arguments de la subroutine rdgrads
	115	c---------------------------------------------------------------------
	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
	130	c
	131	c---------------------------------------------------------------------
	132	c variable argument de la subroutine copie
	133	c---------------------------------------------------------------------
	134	c SB real pplay(100) !pression en milieu de couche du gcm
	135	c SB !argument de la physique
	136	c---------------------------------------------------------------------
	137	c variables destinees a la lecture du pas de temps du fichier de donnees
	138	c---------------------------------------------------------------------
	139	character*80 aaa,atemps,spaces,apasmax
	140	integer nch,imn,ipa
	141	c---------------------------------------------------------------------
	142	c procedures appelees
	143	external rdgrads !lire en iterant dans forcing.dat
	144	c---------------------------------------------------------------------
	145	print*,'le pas itap est:',itap
	146	c* on determine le pas du meso_NH correspondant au nouvel itap *
	147	c* pour aller chercher les champs dans rdgrads *
	148	c
	149	time=time0+itap*dtime
	150	cc temps=int(time/dt+1)
	151	cc pas=min(temps,pasmax)
	152	temps = 1 + int((dt + 2time)/(2dt))
	153	pas=min(temps,pasmax-1)
	154	print*,'le pas Meso est:',pas
	155	c
	156	c
	157	c===================================================================
	158	c
	159	c* on remplit les champs before avec les champs after du pas *
	160	c* precedent en format gcm *
	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
[1960]	171	tsbef = tsaft
[1607]	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
	182	c* on va chercher les nouveaux champs after dans toga.dat *
	183	c* 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	c
	190
	191	if(Tp_fcg) then
	192	c (le forcage est donne en temperature potentielle)
	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
	202	c
	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
[1960]	214	IF (ts_fcg) print*,'ts_subr', ts_subr
[1607]	215	c* on interpole les champs meso_NH sur les niveaux de pression*
	216	c* gcm . on obtient le nouveau champ after *
	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
	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)
	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
	252	c* si on atteint le pas max des donnees experimentales ,on *
	253	c* on conserve les derniers champs calcules *
	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
[1960]	264	ts_subr = tsaft
[1607]	265	else ! temps.ge.pasmax
	266	c* on interpole sur les pas de temps de 10mn du gcm a partir *
	267	c des pas de temps de 1h du meso_NH *
	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
	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
	281	endif ! Turb_fcg
	282	enddo
[1960]	283	ts_subr = ((timeaft-time)tsbef + (time-timebef)tsaft)/dt
[1607]	284	endif ! temps.ge.pasmax
	285	c
	286	print *,' time,timebef,timeaft',time,timebef,timeaft
	287	print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
	288	do j= 1,klev
	289	print *, j,ht(j),htbef(j),htaft(j),
	290	$ hthturb(j),hthturbbef(j),hthturbaft(j)
	291	enddo
	292	print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
	293	do j= 1,klev
	294	print *, j,hq(j),hqbef(j),hqaft(j),
	295	$ hqturb(j),hqturbbef(j),hqturbaft(j)
	296	enddo
	297	c
	298	c-------------------------------------------------------------------
	299	c
	300	IF (Ts_fcg) Ts = Ts_subr
	301	return
	302	c
	303	c-----------------------------------------------------------------------
	304	c on sort les champs de "convergence" pour l instant initial 'in'
	305	c ceci se passe au pas temps itap=0 de la physique
	306	c-----------------------------------------------------------------------
	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)
	310	print*,'le pas itap est:',itap
	311	c
	312	c===================================================================
	313	c
	314	write(*,'(a)') 'OPEN '//file_forctl
	315	open(97,FILE=file_forctl,FORM='FORMATTED')
	316	c
	317	c------------------
	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)
	324	c-------------------------------------------------------------------
	325	c * on lit le pas de temps dans le fichier de donnees *
	326	c * "forcing.ctl" et pasmax *
	327	c-------------------------------------------------------------------
	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)
	346	c------------------------------------------------------------------
	347	c * on lit le pas de temps initial de la simulation *
	348	c------------------------------------------------------------------
	349	in=itap
	350	cc pasprev=in
	351	cc time0=dt*(pasprev-1)
	352	pasprev=in-1
	353	time0=dt*pasprev
	354	C
	355	close(98)
	356	c
	357	write(*,'(a)') 'OPEN '//file_fordat
	358	open(99,FILE=file_fordat,FORM='UNFORMATTED',
	359	. ACCESS='DIRECT',RECL=8)
	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)
	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)
	368	print *, 'get_uvd : rdgrads ->'
	369	print *, tp_fcg
	370	c
	371	c following commented out because we have temperature already in ARM case
	372	c (otherwise this is the potential temperature )
	373	c------------------------------------------------------------------------
	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
	391	c----------------------------------------------------------------------
	392	c on a obtenu des champs initiaux sur les niveaux du meso_NH
	393	c on interpole sur les niveaux du gcm(niveau pression bien sur!)
	394	c-----------------------------------------------------------------------
	395	do k=1,klev
	396	if (JM(k) .eq. 0) then
	397	cFKC bug? ne faut il pas convertir tsol en tendance ????
	398	cRT bug taken care of by removing the stuff
	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
	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)
	423	endif ! Turb_fcg
	424	endif ! JM(k) .eq. 0
	425	enddo
	426	tsaft = ts_subr
	427	c valeurs initiales des champs de convergence
	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
[1960]	441	ts_subr = tsaft
[1607]	442	close(99)
	443	close(98)
	444	c
	445	c-------------------------------------------------------------------
	446	c
	447	c
	448	100 IF (Ts_fcg) Ts = Ts_subr
	449	return
	450	c
	451	999 continue
	452	stop 'erreur lecture, file forcing.ctl'
	453	end
	454
	455	SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f
	456	: ,d_t_adv,d_q_adv)
	457	use dimphy
	458	implicit none
	459
	460	#include "dimensions.h"
	461	ccccc#include "dimphy.h"
	462
	463	integer k
	464	real dtime, fact, du, dv, cx, cy, alx, aly
	465	real zt(klev), zq(klev,3)
	466	: , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
	467
	468	real d_t_adv(klev), d_q_adv(klev,3)
	469
	470	c Velocity of moving cell
	471	data cx,cy /12., -2./
	472
	473	c Dimensions of moving cell
	474	data alx,aly /100 000.,150 000./
	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
	492	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	493	c cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
	494	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	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
[1960]	506	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	507	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
[1607]	508
[1960]	509	integer k,klevgcm
[1607]	510	real playgcm(klevgcm) ! pression en milieu de couche du gcm
	511	real psolgcm
	512	character*80 file_forctl
	513
	514	klev = klevgcm
	515
	516	c---------------------------------------------------------------------
	517	c pression au milieu des couches du gcm dans la physiq
	518	c (SB: remplace le call conv_lipress_gcm(playgcm) )
	519	c---------------------------------------------------------------------
	520
	521	do k = 1, klev
	522	play(k) = playgcm(k)
	523	print*,'la pression gcm est:',play(k)
	524	enddo
	525
	526	c----------------------------------------------------------------------
	527	c lecture du descripteur des donnees Meso-NH (forcing.ctl):
	528	c -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
	529	c (on remplit le COMMON com2_phys_gcss)
	530	c----------------------------------------------------------------------
	531
	532	call mesolupbis(file_forctl)
	533
	534	print*,'la valeur de nblvlm est:',nblvlm
	535
	536	c----------------------------------------------------------------------
	537	c etude de la correspondance entre les niveaux meso.NH et GCM;
	538	c calcul des coefficients d interpolation coef1 et coef2
	539	c (on remplit le COMMON com1_phys_gcss)
	540	c----------------------------------------------------------------------
	541
	542	call corresbis(psolgcm)
	543
	544	c---------------------------------------------------------
	545	c TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
	546	c---------------------------------------------------------
	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
	564	c
	565	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	566	c
	567	c Lecture descripteur des donnees MESO-NH (forcing.ctl):
	568	c -------------------------------------------------------
	569	c
	570	c Cette subroutine lit dans le fichier de controle "essai.ctl"
	571	c et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
	572	c des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
	573	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	574	c
	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
[1960]	578	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
[1607]	579
[1960]	580	INTEGER i,lu,mlz,mlzh
[1607]	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')
	590	c
	591	do i=1,1000
	592	read(lu,1000,end=999) a
	593	if (a .eq. 'ZDEF') go to 100
	594	enddo
	595	c
	596	100 backspace(lu)
	597	print*,' DESCRIPTION DES 2 MODELES : '
	598	print*,' '
	599	c
	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
	606	c
	607	print*,'nbre de niveaux de pression Meso-NH :',nblvlm
	608	print*,' '
	609	print*,'pression en Pa de chaque couche du meso-NH :'
	610	c
	611	read(lu,*) (playm(mlz),mlz=1,nblvlm)
	612	c Si la pression est en HPa, la multiplier par 100
	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)
	619	c
	620	1000 format (a4)
	621	1001 format(5x,i2)
	622	c
	623	print*,' '
	624	do mlzh=1,nblvlm
	625	hplaym(mlzh)=playm(mlzh)/100.
	626	enddo
	627	c
	628	print*,'pression en hPa de chaque couche du meso-NH: '
	629	print*,(hplaym(mlzh),mlzh=1,nblvlm)
	630	c
	631	close (lu)
	632	return
	633	c
	634	999 stop 'erreur lecture des niveaux pression des donnees'
	635	end
	636
	637	SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,
	638	: ts_fcg,ts,imp_fcg,Turb_fcg)
	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
	644	c
[1960]	645	INTEGER k
[1607]	646	c
	647	LOGICAL imp_fcg,ts_fcg,Turb_fcg
	648	c
	649	icomp = icount
	650	c
	651	c
	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
	666	c
	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'
	678	c
	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
	691	c
	692	do k=1,nl
	693	icomp=icomp+1
	694	read(itape,rec=icomp)hw(k)
	695	enddo
	696	c
	697	if(ts_fcg) then
	698	icomp=icomp+1
	699	read(itape,rec=icomp)ts
	700	endif
	701	c
	702	print *,' rdgrads ->'
	703
	704	RETURN
	705	END
	706
	707	SUBROUTINE corresbis(psol)
	708	implicit none
	709
	710	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	711	c Cette subroutine calcule et affiche les valeurs des coefficients
	712	c d interpolation qui serviront dans la formule d interpolation elle-
	713	c meme.
	714	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	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
[1960]	726	COMMON/com1_phys_gcss/play,coef1,coef2,JM,klev
	727	COMMON/com2_phys_gcss/playm,hplaym,nblvlm
[1607]	728
	729	REAL psol
	730	REAL val
[1960]	731	INTEGER k, mlz
[1607]	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
	739	coef1(1)=(playm(mlz+1)-val)
	740	* /(playm(mlz+1)-psol)
	741	coef2(1)=(val-psol)
	742	* /(playm(mlz+1)-psol)
	743	else if (val .gt. playm(nblvlm)) then
	744	do mlz=1,nblvlm
	745	if ( val .le. playm(mlz)
	746	* .and. val .gt. playm(mlz+1))then
	747	JM(k)=mlz
	748	coef1(k)=(playm(mlz+1)-val)
	749	* /(playm(mlz+1)-playm(mlz))
	750	coef2(k)=(val-playm(mlz))
	751	* /(playm(mlz+1)-playm(mlz))
	752	endif
	753	enddo
	754	else
	755	JM(k) = nblvlm-1
	756	coef1(k) = 0.
	757	coef2(k) = 0.
	758	endif
	759	enddo
	760	c
	761	cc if (play(klev) .le. playm(nblvlm)) then
	762	cc mlz=nblvlm-1
	763	cc JM(klev)=mlz
	764	cc coef1(klev)=(playm(mlz+1)-val)
	765	cc * /(playm(mlz+1)-playm(mlz))
	766	cc coef2(klev)=(val-playm(mlz))
	767	cc * /(playm(mlz+1)-playm(mlz))
	768	cc endif
	769	c
	770	print*,' '
	771	print*,' INTERPOLATION : '
	772	print*,' '
	773	print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
	774	print*,(JM(k),k=1,klev)
	775	print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
	776	print*,(JM(k),k=1,klev)
	777	print*,' '
	778	print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
	779	*: '
	780	print*,(coef1(k),k=1,klev)
	781	print*,' '
	782	print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
	783	*x: '
	784	print*,(coef2(k),k=1,klev)
	785	c
	786	return
	787	end
	788	SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
	789	C***************************************************************
	790	C* *
	791	C* *
	792	C* GETSCH *
	793	C* *
	794	C* *
	795	C* modified by : *
	796	C***************************************************************
	797	C* Return in SST the character string found between the NTH-1 and NTH
	798	C* occurence of the delimiter 'DEL' but before the terminator 'TRM' in
	799	C* the input string 'STR'. If TRM=DEL then STR is considered unlimited.
	800	C* NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
	801	C* NTH is greater than the number of delimiters in STR.
	802	IMPLICIT INTEGER (A-Z)
	803	CHARACTER STR(),DEL1,TRM1,SST()
	804	NCH=-1
	805	SST=' '
	806	IF(NTH.GT.0) THEN
	807	IF(TRM.EQ.DEL) THEN
	808	LENGTH=LEN(STR)
	809	ELSE
	810	LENGTH=INDEX(STR,TRM)-1
	811	IF(LENGTH.LT.0) LENGTH=LEN(STR)
	812	ENDIF
	813	C* Find beginning and end of the NTH DEL-limited substring in STR
	814	END=-1
	815	DO 1,N=1,NTH
	816	IF(END.EQ.LENGTH) RETURN
	817	BEG=END+2
	818	END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
	819	IF(END.EQ.BEG-2) END=LENGTH
	820	C* PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
	821	1 CONTINUE
	822	NCH=END-BEG+1
	823	IF(NCH.GT.0) SST=STR(BEG:END)
	824	ENDIF
	825	END
	826	CHARACTER() FUNCTION SPACES(STR,NSPACE)
	827	C
	828	C CERN PROGLIB# M433 SPACES .VERSION KERNFOR 4.14 860211
	829	C ORIG. 6/05/86 M.GOOSSENS/DD
	830	C
	831	C- The function value SPACES returns the character string STR with
	832	C- leading blanks removed and each occurence of one or more blanks
	833	C- replaced by NSPACE blanks inside the string STR
	834	C
	835	CHARACTER() STR
	836	C
	837	LENSPA = LEN(SPACES)
	838	SPACES = ' '
	839	IF (NSPACE.LT.0) NSPACE = 0
	840	IBLANK = 1
	841	ISPACE = 1
	842	100 INONBL = INDEXC(STR(IBLANK:),' ')
	843	IF (INONBL.EQ.0) THEN
	844	SPACES(ISPACE:) = STR(IBLANK:)
	845	GO TO 999
	846	ENDIF
	847	INONBL = INONBL + IBLANK - 1
	848	IBLANK = INDEX(STR(INONBL:),' ')
	849	IF (IBLANK.EQ.0) THEN
	850	SPACES(ISPACE:) = STR(INONBL:)
	851	GO TO 999
	852	ENDIF
	853	IBLANK = IBLANK + INONBL - 1
	854	SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
	855	ISPACE = ISPACE + IBLANK - INONBL + NSPACE
	856	IF (ISPACE.LE.LENSPA) GO TO 100
	857	999 END
	858	FUNCTION INDEXC(STR,SSTR)
	859	C
	860	C CERN PROGLIB# M433 INDEXC .VERSION KERNFOR 4.14 860211
	861	C ORIG. 26/03/86 M.GOOSSENS/DD
	862	C
	863	C- Find the leftmost position where substring SSTR does not match
	864	C- string STR scanning forward
	865	C
	866	CHARACTER() STR,SSTR
	867	C
	868	LENS = LEN(STR)
	869	LENSS = LEN(SSTR)
	870	C
	871	DO 10 I=1,LENS-LENSS+1
	872	IF (STR(I:I+LENSS-1).NE.SSTR) THEN
	873	INDEXC = I
	874	GO TO 999
	875	ENDIF
	876	10 CONTINUE
	877	INDEXC = 0
	878	C
	879	999 END

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

Download in other formats: