Context Navigation

calfis.F @ 5306

Last change on this file since 5306 was 344, checked in by lmdz, 23 years ago
Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 16.6 KB

Rev	Line
[344]	1	SUBROUTINE calfis(nq,
	2	$ lafin,
	3	$ rdayvrai,
	4	$ rday_ecri,
	5	$ heure,
	6	$ pucov,
	7	$ pvcov,
	8	$ pteta,
	9	$ pq,
	10	$ pmasse,
	11	$ pps,
	12	$ pp,
	13	$ ppk,
	14	$ pphis,
	15	$ pphi,
	16	$ pducov,
	17	$ pdvcov,
	18	$ pdteta,
	19	$ pdq,
	20	$ pw,
	21	#ifdef INCA_CH4
	22	$ flxw,
	23	#endif
	24	$ clesphy0,
	25	$ pdufi,
	26	$ pdvfi,
	27	$ pdhfi,
	28	$ pdqfi,
	29	$ pdpsfi)
[2]	30	c
	31	c Auteur : P. Le Van, F. Hourdin
	32	c .........
	33
	34	IMPLICIT NONE
	35	c=======================================================================
	36	c
	37	c 1. rearrangement des tableaux et transformation
	38	c variables dynamiques > variables physiques
	39	c 2. calcul des termes physiques
	40	c 3. retransformation des tendances physiques en tendances dynamiques
	41	c
	42	c remarques:
	43	c ----------
	44	c
	45	c - les vents sont donnes dans la physique par leurs composantes
	46	c naturelles.
	47	c - la variable thermodynamique de la physique est une variable
	48	c intensive : T
	49	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
	50	c - les deux seules variables dependant de la geometrie necessaires
	51	c pour la physique sont la latitude pour le rayonnement et
	52	c l'aire de la maille quand on veut integrer une grandeur
	53	c horizontalement.
	54	c - les points de la physique sont les points scalaires de la
	55	c la dynamique; numerotation:
	56	c 1 pour le pole nord
	57	c (jjm-1)*iim pour l'interieur du domaine
	58	c ngridmx pour le pole sud
	59	c ---> ngridmx=2+(jjm-1)*iim
	60	c
	61	c Input :
	62	c -------
	63	c ecritphy frequence d'ecriture (en jours)de histphy
	64	c pucov covariant zonal velocity
	65	c pvcov covariant meridional velocity
	66	c pteta potential temperature
	67	c pps surface pressure
	68	c pmasse masse d'air dans chaque maille
	69	c pts surface temperature (K)
	70	c callrad clef d'appel au rayonnement
	71	c
	72	c Output :
	73	c --------
	74	c pdufi tendency for the natural zonal velocity (ms-1)
	75	c pdvfi tendency for the natural meridional velocity
	76	c pdhfi tendency for the potential temperature
	77	c pdtsfi tendency for the surface temperature
	78	c
	79	c pdtrad radiative tendencies \ both input
	80	c pfluxrad radiative fluxes / and output
	81	c
	82	c=======================================================================
	83	c
	84	c-----------------------------------------------------------------------
	85	c
	86	c 0. Declarations :
	87	c ------------------
	88
	89	#include "dimensions.h"
	90	#include "paramet.h"
	91	#include "temps.h"
	92
	93	INTEGER ngridmx,nq
	94	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	95
	96	#include "comconst.h"
	97	#include "comvert.h"
	98	#include "comgeom2.h"
	99	#include "control.h"
	100
	101	c Arguments :
	102	c -----------
	103	LOGICAL lafin
	104	REAL heure
	105
	106	REAL pvcov(iip1,jjm,llm)
	107	REAL pucov(iip1,jjp1,llm)
	108	REAL pteta(iip1,jjp1,llm)
	109	REAL pmasse(iip1,jjp1,llm)
	110	REAL pq(iip1,jjp1,llm,nqmx)
	111	REAL pphis(iip1,jjp1)
	112	REAL pphi(iip1,jjp1,llm)
	113	c
	114	REAL pdvcov(iip1,jjm,llm)
	115	REAL pducov(iip1,jjp1,llm)
	116	REAL pdteta(iip1,jjp1,llm)
	117	REAL pdq(iip1,jjp1,llm,nqmx)
	118	c
	119	REAL pw(iip1,jjp1,llm)
	120	c
	121	REAL pps(iip1,jjp1)
	122	REAL pp(iip1,jjp1,llmp1)
	123	REAL ppk(iip1,jjp1,llm)
	124	c
	125	REAL pdvfi(iip1,jjm,llm)
	126	REAL pdufi(iip1,jjp1,llm)
	127	REAL pdhfi(iip1,jjp1,llm)
	128	REAL pdqfi(iip1,jjp1,llm,nqmx)
	129	REAL pdpsfi(iip1,jjp1)
	130
	131	INTEGER longcles
	132	PARAMETER ( longcles = 20 )
	133	REAL clesphy0( longcles )
	134
	135
	136	c Local variables :
	137	c -----------------
	138
	139	INTEGER i,j,l,ig0,ig,iq
	140	REAL zpsrf(ngridmx)
	141	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
	142	REAL zphi(ngridmx,llm),zphis(ngridmx)
	143	c
	144	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
	145	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
	146	c
	147	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
	148	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
	149	c
	150	REAL pvervel(ngridmx,llm)
	151	c
	152	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
	153	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
	154	REAL zdpsrf(ngridmx)
	155	c
	156	REAL zsin(iim),zcos(iim),z1(iim)
	157	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	158	REAL unskap, pksurcp
[344]	159
	160	#ifdef INCA_CH4
	161	REAL flxw(iip1,jjp1,llm)
	162	REAL flxwfi(ngridmx,llm)
	163	#endif
[2]	164
	165	EXTERNAL gr_dyn_fi,gr_fi_dyn
	166	EXTERNAL physiq,multipl
	167	REAL SSUM
	168	EXTERNAL SSUM
	169
	170	REAL latfi(ngridmx),lonfi(ngridmx)
	171	REAL airefi(ngridmx)
[18]	172	SAVE latfi, lonfi, airefi
[2]	173
	174	LOGICAL firstcal, debut
	175	DATA firstcal/.true./
	176	SAVE firstcal,debut
	177	REAL rdayvrai,rday_ecri
	178	c
	179	c-----------------------------------------------------------------------
	180	c
	181	c 1. Initialisations :
	182	c --------------------
	183	c
	184
	185	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
	186	PRINT*,'STOP dans calfis'
	187	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
	188	PRINT*,' ngridmx jjm iim '
	189	PRINT*,ngridmx,jjm,iim
	190	STOP
	191	ENDIF
	192
	193	c-----------------------------------------------------------------------
	194	c latitude, longitude et aires des mailles pour la physique:
	195	c ----------------------------------------------------------
	196
	197	c
	198	IF ( firstcal ) THEN
	199	debut = .TRUE.
	200	ELSE
	201	debut = .FALSE.
	202	ENDIF
	203
	204	c
	205	IF (firstcal) THEN
	206	latfi(1)=rlatu(1)
	207	lonfi(1)=0.
	208	DO j=2,jjm
	209	DO i=1,iim
	210	latfi((j-2)*iim+1+i)= rlatu(j)
	211	lonfi((j-2)*iim+1+i)= rlonv(i)
	212	ENDDO
	213	ENDDO
	214	latfi(ngridmx)= rlatu(jjp1)
	215	lonfi(ngridmx)= 0.
	216	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
	217	PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
	218	CALL inifis(ngridmx,llm,daysec,day_ini,dtphys ,
	219	, latfi,lonfi,airefi,rad,g,r,cpp )
	220	ENDIF
	221
	222	c
	223	c-----------------------------------------------------------------------
	224	c 40. transformation des variables dynamiques en variables physiques:
	225	c ---------------------------------------------------------------
	226
	227	c 41. pressions au sol (en Pascals)
	228	c ----------------------------------
	229
	230
	231	zpsrf(1) = pps(1,1)
	232
	233	ig0 = 2
	234	DO j = 2,jjm
	235	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
	236	ig0 = ig0+iim
	237	ENDDO
	238
	239	zpsrf(ngridmx) = pps(1,jjp1)
	240
	241
	242	c 42. pression intercouches :
	243	c
	244	c -----------------------------------------------------------------
	245	c .... zplev definis aux (llm +1) interfaces des couches ....
	246	c .... zplay definis aux ( llm ) milieux des couches ....
	247	c -----------------------------------------------------------------
	248
	249	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
	250	c
	251	unskap = 1./ kappa
	252	c
	253	DO l = 1, llmp1
	254	zplev( 1,l ) = pp(1,1,l)
	255	ig0 = 2
	256	DO j = 2, jjm
	257	DO i =1, iim
	258	zplev( ig0,l ) = pp(i,j,l)
	259	ig0 = ig0 +1
	260	ENDDO
	261	ENDDO
	262	zplev( ngridmx,l ) = pp(1,jjp1,l)
	263	ENDDO
	264	c
	265	c
	266
	267	c 43. temperature naturelle (en K) et pressions milieux couches .
	268	c ---------------------------------------------------------------
	269
	270	DO l=1,llm
	271
	272	pksurcp = ppk(1,1,l) / cpp
	273	zplay(1,l) = preff * pksurcp ** unskap
	274	ztfi(1,l) = pteta(1,1,l) * pksurcp
	275	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
	276	ig0 = 2
	277
	278	DO j = 2, jjm
	279	DO i = 1, iim
	280	pksurcp = ppk(i,j,l) / cpp
	281	zplay(ig0,l) = preff * pksurcp ** unskap
	282	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
	283	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
	284	ig0 = ig0 + 1
	285	ENDDO
	286	ENDDO
	287
	288	pksurcp = ppk(1,jjp1,l) / cpp
	289	zplay(ig0,l) = preff * pksurcp ** unskap
	290	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
	291	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
	292
	293	ENDDO
	294
	295	c 43.bis humidite specifique (en kg/kg)
	296	c -------------------------------------
	297
	298	DO iq=1,nqmx
	299	DO l=1,llm
	300	zqfi(1,l,iq) = pq(1,1,l,iq)
	301	ig0 = 2
	302	DO j=2,jjm
	303	DO i = 1, iim
	304	zqfi(ig0,l,iq) = pq(i,j,l,iq)
	305	ig0 = ig0 + 1
	306	ENDDO
	307	ENDDO
	308	zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
	309	ENDDO
	310	ENDDO
	311
	312	c convergence dynamique pour les traceurs "EAU"
	313
	314	DO iq=1,2
	315	DO l=1,llm
	316	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
	317	ig0 = 2
	318	DO j=2,jjm
	319	DO i = 1, iim
	320	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
	321	ig0 = ig0 + 1
	322	ENDDO
	323	ENDDO
	324	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
	325	ENDDO
	326	ENDDO
	327
	328
	329	c Geopotentiel calcule par rapport a la surface locale:
	330	c -----------------------------------------------------
	331
	332	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
	333	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
	334	DO l=1,llm
	335	DO ig=1,ngridmx
	336	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	337	ENDDO
	338	ENDDO
	339
	340	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
	341	c
	342	DO l=1,llm
[343]	343	pvervel(1,l)=pw(1,1,l) /apoln
[2]	344	ig0=2
	345	DO j=2,jjm
	346	DO i = 1, iim
[343]	347	pvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
[2]	348	ig0 = ig0 + 1
	349	ENDDO
	350	ENDDO
[343]	351	pvervel(ig0,l)=pw(1,jjp1,l) /apols
[2]	352	ENDDO
	353
	354	c CALL initial0( ngridmx*llm,pvervel )
	355	c
	356	c
	357	c 45. champ u:
	358	c ------------
	359
	360	DO 50 l=1,llm
	361
	362	DO 25 j=2,jjm
	363	ig0 = 1+(j-2)*iim
	364	zufi(ig0+1,l)= 0.5 *
	365	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
	366	pcvgu(ig0+1,l)= 0.5 *
	367	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
	368	DO 10 i=2,iim
	369	zufi(ig0+i,l)= 0.5 *
	370	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
	371	pcvgu(ig0+i,l)= 0.5 *
	372	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
	373	10 CONTINUE
	374	25 CONTINUE
	375
	376	50 CONTINUE
	377
	378
	379	c 46.champ v:
	380	c -----------
	381
	382	DO l=1,llm
	383	DO j=2,jjm
	384	ig0=1+(j-2)*iim
	385	DO i=1,iim
	386	zvfi(ig0+i,l)= 0.5 *
	387	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
	388	pcvgv(ig0+i,l)= 0.5 *
	389	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
	390	ENDDO
	391	ENDDO
	392	ENDDO
	393
	394
	395	c 47. champs de vents aux pole nord
	396	c ------------------------------
	397	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	398	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	399
	400	DO l=1,llm
	401
	402	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	403	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
	404	DO i=2,iim
	405	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	406	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
	407	ENDDO
	408
	409	DO i=1,iim
	410	zcos(i) = COS(rlonv(i))*z1(i)
	411	zcosbis(i)= COS(rlonv(i))*z1bis(i)
	412	zsin(i) = SIN(rlonv(i))*z1(i)
	413	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
	414	ENDDO
	415
	416	zufi(1,l) = SSUM(iim,zcos,1)/pi
	417	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
	418	zvfi(1,l) = SSUM(iim,zsin,1)/pi
	419	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
	420
	421	ENDDO
	422
	423
	424	c 48. champs de vents aux pole sud:
	425	c ---------------------------------
	426	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	427	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	428
	429	DO l=1,llm
	430
	431	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	432	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
	433	DO i=2,iim
	434	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	435	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
	436	ENDDO
	437
	438	DO i=1,iim
	439	zcos(i) = COS(rlonv(i))*z1(i)
	440	zcosbis(i) = COS(rlonv(i))*z1bis(i)
	441	zsin(i) = SIN(rlonv(i))*z1(i)
	442	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
	443	ENDDO
	444
	445	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
	446	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
	447	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
	448	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
	449
	450	ENDDO
	451
[344]	452	#ifdef INCA_CH4
	453	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
	454	#endif
	455
[2]	456	c-----------------------------------------------------------------------
	457	c Appel de la physique:
	458	c ---------------------
	459
	460
[344]	461	CALL physiq (ngridmx,
	462	. llm,
	463	. nq,
	464	. debut,
	465	. lafin,
	466	. rdayvrai,
	467	. rday_ecri,
	468	. heure,
	469	. dtphys,
	470	. zplev,
	471	. zplay,
	472	. zphi,
	473	. zphis,
	474	. airefi,
	475	. presnivs,
	476	. clesphy0,
	477	. zufi,
	478	. zvfi,
	479	. ztfi,
	480	. zqfi,
	481	. pvervel,
	482	#ifdef INCA_CH4
	483	. flxwfi,
	484	#endif
	485	. zdufi,
	486	. zdvfi,
	487	. zdtfi,
	488	. zdqfi,
	489	. zdpsrf)
[2]	490
	491	500 CONTINUE
	492
	493	c-----------------------------------------------------------------------
	494	c transformation des tendances physiques en tendances dynamiques:
	495	c ---------------------------------------------------------------
	496
	497	c tendance sur la pression :
	498	c -----------------------------------
	499
	500	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
	501	c
	502	ccc CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
	503
	504	c 62. enthalpie potentielle
	505	c ---------------------
	506
	507	DO l=1,llm
	508
	509	DO i=1,iip1
	510	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	511	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
	512	ENDDO
	513
	514	DO j=2,jjm
	515	ig0=1+(j-2)*iim
	516	DO i=1,iim
	517	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
	518	ENDDO
	519	pdhfi(iip1,j,l) = pdhfi(1,j,l)
	520	ENDDO
	521
	522	ENDDO
	523
	524
	525	c 62. humidite specifique
	526	c ---------------------
	527
	528	DO iq=1,nqmx
	529	DO l=1,llm
	530	DO i=1,iip1
	531	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	532	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
	533	ENDDO
	534	DO j=2,jjm
	535	ig0=1+(j-2)*iim
	536	DO i=1,iim
	537	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
	538	ENDDO
	539	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
	540	ENDDO
	541	ENDDO
	542	ENDDO
	543
	544	c 65. champ u:
	545	c ------------
	546
	547	DO l=1,llm
	548
	549	DO i=1,iip1
	550	pdufi(i,1,l) = 0.
	551	pdufi(i,jjp1,l) = 0.
	552	ENDDO
	553
	554	DO j=2,jjm
	555	ig0=1+(j-2)*iim
	556	DO i=1,iim-1
	557	pdufi(i,j,l)=
	558	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
	559	ENDDO
	560	pdufi(iim,j,l)=
	561	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
	562	pdufi(iip1,j,l)=pdufi(1,j,l)
	563	ENDDO
	564
	565	ENDDO
	566
	567
	568	c 67. champ v:
	569	c ------------
	570
	571	DO l=1,llm
	572
	573	DO j=2,jjm-1
	574	ig0=1+(j-2)*iim
	575	DO i=1,iim
	576	pdvfi(i,j,l)=
	577	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
	578	ENDDO
	579	pdvfi(iip1,j,l) = pdvfi(1,j,l)
	580	ENDDO
	581	ENDDO
	582
	583
	584	c 68. champ v pres des poles:
	585	c ---------------------------
	586	c v = U * cos(long) + V * SIN(long)
	587
	588	DO l=1,llm
	589
	590	DO i=1,iim
	591	pdvfi(i,1,l)=
	592	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	593	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
	594	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
	595	pdvfi(i,1,l)=
	596	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	597	pdvfi(i,jjm,l)=
	598	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
	599	ENDDO
	600
	601	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	602	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	603
	604	ENDDO
	605
	606	c-----------------------------------------------------------------------
	607
	608	700 CONTINUE
	609
	610	firstcal = .FALSE.
	611
	612	RETURN
	613	END

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

Download in other formats:

Original Format