Context Navigation

calfis.F @ 729

Last change on this file since 729 was 644, checked in by Laurent Fairhead, 19 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.6 KB

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

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

Download in other formats:

Original Format