Context Navigation

calfis.F @ 2192

Last change on this file since 2192 was 2037, checked in by lguez, 11 years ago

PVteta computed by PVtheta was not used. Also there were a couple of
problems with PVtheta:

-- PVtheta calls tetalevel in phylmd, and interpolates at

Earth-specific values of potential temperature.

-- PVtheta calls tourabs, and the computation of rot in tourabs should

be modified (it is not correct when there is a zoom).

-- Even when there is no zoom, the computation of rot in tourabs

should probably be modified: directly combine vcov and ucov, as in
tourpot, instead of dividing by cv and cu.

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:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 19.9 KB

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

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

Download in other formats:

Original Format