Context Navigation

calfis.F @ 1541

Last change on this file since 1541 was 1407, checked in by Laurent Fairhead, 14 years ago

Some extraneous prints deleted for optimisation purposes

Des prints supprimés pour optimisation

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 18.1 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
[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 .........
[1146]	32	USE infotrac
[1403]	33	USE control_mod
	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 -----------
	104	LOGICAL lafin
	105
[1279]	106
[524]	107	REAL pvcov(iip1,jjm,llm)
	108	REAL pucov(iip1,jjp1,llm)
	109	REAL pteta(iip1,jjp1,llm)
	110	REAL pmasse(iip1,jjp1,llm)
[1146]	111	REAL pq(iip1,jjp1,llm,nqtot)
[524]	112	REAL pphis(iip1,jjp1)
	113	REAL pphi(iip1,jjp1,llm)
	114	c
	115	REAL pdvcov(iip1,jjm,llm)
	116	REAL pducov(iip1,jjp1,llm)
	117	REAL pdteta(iip1,jjp1,llm)
[1146]	118	REAL pdq(iip1,jjp1,llm,nqtot)
[524]	119	c
	120	REAL pps(iip1,jjp1)
	121	REAL pp(iip1,jjp1,llmp1)
	122	REAL ppk(iip1,jjp1,llm)
	123	c
	124	REAL pdvfi(iip1,jjm,llm)
	125	REAL pdufi(iip1,jjp1,llm)
	126	REAL pdhfi(iip1,jjp1,llm)
[1146]	127	REAL pdqfi(iip1,jjp1,llm,nqtot)
[524]	128	REAL pdpsfi(iip1,jjp1)
	129
	130	INTEGER longcles
	131	PARAMETER ( longcles = 20 )
	132	REAL clesphy0( longcles )
	133
	134
	135	c Local variables :
	136	c -----------------
	137
	138	INTEGER i,j,l,ig0,ig,iq,iiq
	139	REAL zpsrf(ngridmx)
	140	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
	141	REAL zphi(ngridmx,llm),zphis(ngridmx)
	142	c
	143	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
[1146]	144	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
[524]	145	c
	146	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
	147	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
	148	c
	149	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
[1146]	150	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
[524]	151	REAL zdpsrf(ngridmx)
	152	c
[1403]	153	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
	154	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
	155	REAL jH_cur_split,zdt_split
	156	LOGICAL debut_split,lafin_split
	157	INTEGER isplit
	158
[524]	159	REAL zsin(iim),zcos(iim),z1(iim)
	160	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	161	REAL unskap, pksurcp
[644]	162	c
	163	cIM diagnostique PVteta, Amip2
	164	INTEGER ntetaSTD
	165	PARAMETER(ntetaSTD=3)
	166	REAL rtetaSTD(ntetaSTD)
	167	DATA rtetaSTD/350., 380., 405./
	168	REAL PVteta(ngridmx,ntetaSTD)
	169	c
[960]	170	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
	171	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
[524]	172	c
	173
	174	REAL SSUM
	175
	176	LOGICAL firstcal, debut
	177	DATA firstcal/.true./
	178	SAVE firstcal,debut
[1279]	179	! REAL rdayvrai
	180	REAL, intent(in):: jD_cur, jH_cur
[524]	181	c
	182	c-----------------------------------------------------------------------
	183	c
	184	c 1. Initialisations :
	185	c --------------------
	186	c
[1279]	187	c
	188	IF ( firstcal ) THEN
	189	debut = .TRUE.
	190	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
[1403]	191	write(lunout,*) 'STOP dans calfis'
	192	write(lunout,*)
	193	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
	194	write(lunout,*) ' ngridmx jjm iim '
	195	write(lunout,*) ngridmx,jjm,iim
[524]	196	STOP
[1279]	197	ENDIF
[524]	198	ELSE
[1279]	199	debut = .FALSE.
	200	ENDIF ! of IF (firstcal)
[524]	201
	202	c
	203	c
	204	c-----------------------------------------------------------------------
	205	c 40. transformation des variables dynamiques en variables physiques:
	206	c ---------------------------------------------------------------
	207
	208	c 41. pressions au sol (en Pascals)
	209	c ----------------------------------
	210
	211
	212	zpsrf(1) = pps(1,1)
	213
	214	ig0 = 2
	215	DO j = 2,jjm
	216	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
	217	ig0 = ig0+iim
	218	ENDDO
	219
	220	zpsrf(ngridmx) = pps(1,jjp1)
	221
	222
	223	c 42. pression intercouches :
	224	c
	225	c -----------------------------------------------------------------
	226	c .... zplev definis aux (llm +1) interfaces des couches ....
	227	c .... zplay definis aux ( llm ) milieux des couches ....
	228	c -----------------------------------------------------------------
	229
	230	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
	231	c
	232	unskap = 1./ kappa
	233	c
	234	DO l = 1, llmp1
	235	zplev( 1,l ) = pp(1,1,l)
	236	ig0 = 2
	237	DO j = 2, jjm
	238	DO i =1, iim
	239	zplev( ig0,l ) = pp(i,j,l)
	240	ig0 = ig0 +1
	241	ENDDO
	242	ENDDO
	243	zplev( ngridmx,l ) = pp(1,jjp1,l)
	244	ENDDO
	245	c
	246	c
	247
	248	c 43. temperature naturelle (en K) et pressions milieux couches .
	249	c ---------------------------------------------------------------
	250
	251	DO l=1,llm
	252
	253	pksurcp = ppk(1,1,l) / cpp
	254	zplay(1,l) = preff * pksurcp ** unskap
	255	ztfi(1,l) = pteta(1,1,l) * pksurcp
	256	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
	257	ig0 = 2
	258
	259	DO j = 2, jjm
	260	DO i = 1, iim
	261	pksurcp = ppk(i,j,l) / cpp
	262	zplay(ig0,l) = preff * pksurcp ** unskap
	263	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
	264	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
	265	ig0 = ig0 + 1
	266	ENDDO
	267	ENDDO
	268
	269	pksurcp = ppk(1,jjp1,l) / cpp
	270	zplay(ig0,l) = preff * pksurcp ** unskap
	271	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
	272	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
	273
	274	ENDDO
	275
	276	c 43.bis traceurs
	277	c ---------------
	278	c
[1146]	279	DO iq=1,nqtot
[524]	280	iiq=niadv(iq)
	281	DO l=1,llm
	282	zqfi(1,l,iq) = pq(1,1,l,iiq)
	283	ig0 = 2
	284	DO j=2,jjm
	285	DO i = 1, iim
	286	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
	287	ig0 = ig0 + 1
	288	ENDDO
	289	ENDDO
	290	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
	291	ENDDO
	292	ENDDO
	293
	294	c convergence dynamique pour les traceurs "EAU"
[1279]	295	! Earth-specific treatment of first 2 tracers (water)
	296	if (planet_type=="earth") then
	297	DO iq=1,2
[524]	298	DO l=1,llm
	299	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
	300	ig0 = 2
	301	DO j=2,jjm
	302	DO i = 1, iim
	303	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
	304	ig0 = ig0 + 1
	305	ENDDO
	306	ENDDO
	307	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
	308	ENDDO
[1279]	309	ENDDO
	310	endif ! of if (planet_type=="earth")
[524]	311
	312
	313	c Geopotentiel calcule par rapport a la surface locale:
	314	c -----------------------------------------------------
	315
	316	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
	317	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
	318	DO l=1,llm
[1403]	319	DO ig=1,ngridmx
	320	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	321	ENDDO
[524]	322	ENDDO
	323
	324	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
[960]	325	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
	326	c de masse est calclue dans advtrac.F
	327	c DO l=1,llm
	328	c pvervel(1,l)=pw(1,1,l) * g /apoln
	329	c ig0=2
	330	c DO j=2,jjm
	331	c DO i = 1, iim
	332	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
	333	c ig0 = ig0 + 1
	334	c ENDDO
	335	c ENDDO
	336	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
	337	c ENDDO
[524]	338
	339	c
	340	c 45. champ u:
	341	c ------------
	342
	343	DO 50 l=1,llm
	344
	345	DO 25 j=2,jjm
	346	ig0 = 1+(j-2)*iim
	347	zufi(ig0+1,l)= 0.5 *
	348	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
	349	pcvgu(ig0+1,l)= 0.5 *
	350	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
	351	DO 10 i=2,iim
	352	zufi(ig0+i,l)= 0.5 *
	353	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
	354	pcvgu(ig0+i,l)= 0.5 *
	355	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
	356	10 CONTINUE
	357	25 CONTINUE
	358
	359	50 CONTINUE
	360
	361
	362	c 46.champ v:
	363	c -----------
	364
	365	DO l=1,llm
	366	DO j=2,jjm
	367	ig0=1+(j-2)*iim
	368	DO i=1,iim
	369	zvfi(ig0+i,l)= 0.5 *
	370	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
	371	pcvgv(ig0+i,l)= 0.5 *
	372	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
	373	ENDDO
	374	ENDDO
	375	ENDDO
	376
	377
	378	c 47. champs de vents aux pole nord
	379	c ------------------------------
	380	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	381	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	382
	383	DO l=1,llm
	384
	385	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	386	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
	387	DO i=2,iim
	388	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	389	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
	390	ENDDO
	391
	392	DO i=1,iim
	393	zcos(i) = COS(rlonv(i))*z1(i)
	394	zcosbis(i)= COS(rlonv(i))*z1bis(i)
	395	zsin(i) = SIN(rlonv(i))*z1(i)
	396	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
	397	ENDDO
	398
	399	zufi(1,l) = SSUM(iim,zcos,1)/pi
	400	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
	401	zvfi(1,l) = SSUM(iim,zsin,1)/pi
	402	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
	403
	404	ENDDO
	405
	406
	407	c 48. champs de vents aux pole sud:
	408	c ---------------------------------
	409	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	410	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	411
	412	DO l=1,llm
	413
	414	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	415	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
	416	DO i=2,iim
	417	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	418	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
[1403]	419	ENDDO
[524]	420
	421	DO i=1,iim
	422	zcos(i) = COS(rlonv(i))*z1(i)
	423	zcosbis(i) = COS(rlonv(i))*z1bis(i)
	424	zsin(i) = SIN(rlonv(i))*z1(i)
	425	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
[1403]	426	ENDDO
[524]	427
	428	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
	429	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
	430	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
	431	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
	432
	433	ENDDO
[644]	434	c
[1279]	435	if (planet_type=="earth") then
	436	#ifdef CPP_EARTH
[644]	437	cIM calcul PV a teta=350, 380, 405K
	438	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
	439	$ ztfi,zplay,zplev,
	440	$ ntetaSTD,rtetaSTD,PVteta)
[1279]	441	#endif
	442	endif
[644]	443	c
[960]	444	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
[524]	445	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
	446
	447	c-----------------------------------------------------------------------
	448	c Appel de la physique:
	449	c ---------------------
	450
	451
[1279]	452	if (planet_type=="earth") then
	453	#ifdef CPP_EARTH
[1403]	454
[1407]	455	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
[1403]	456	zdt_split=dtphys/nsplit_phys
	457	zdufic(:,:)=0.
	458	zdvfic(:,:)=0.
	459	zdtfic(:,:)=0.
	460	zdqfic(:,:,:)=0.
	461
	462	do isplit=1,nsplit_phys
	463
	464	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
	465	debut_split=debut.and.isplit==1
	466	lafin_split=lafin.and.isplit==nsplit_phys
	467
	468	CALL physiq (ngridmx,
[524]	469	. llm,
[1403]	470	. debut_split,
	471	. lafin_split,
[1279]	472	. jD_cur,
[1403]	473	. jH_cur_split,
	474	. zdt_split,
[524]	475	. zplev,
	476	. zplay,
	477	. zphi,
	478	. zphis,
	479	. presnivs,
	480	. clesphy0,
	481	. zufi,
	482	. zvfi,
	483	. ztfi,
	484	. zqfi,
	485	. flxwfi,
	486	. zdufi,
	487	. zdvfi,
	488	. zdtfi,
	489	. zdqfi,
[644]	490	. zdpsrf,
	491	cIM diagnostique PVteta, Amip2
	492	. pducov,
	493	. PVteta)
[1403]	494
	495	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
	496	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
	497	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
	498	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
	499
	500	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
	501	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
	502	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
	503	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
	504
	505	enddo
	506	zdufi(:,:)=zdufic(:,:)/nsplit_phys
	507	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
	508	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
	509	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
	510
[1279]	511	#endif
	512	endif !of if (planet_type=="earth")
[524]	513
	514	500 CONTINUE
	515
	516	c-----------------------------------------------------------------------
	517	c transformation des tendances physiques en tendances dynamiques:
	518	c ---------------------------------------------------------------
	519
	520	c tendance sur la pression :
	521	c -----------------------------------
	522
	523	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
	524	c
	525	c 62. enthalpie potentielle
	526	c ---------------------
	527
	528	DO l=1,llm
	529
	530	DO i=1,iip1
	531	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	532	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
	533	ENDDO
	534
	535	DO j=2,jjm
	536	ig0=1+(j-2)*iim
	537	DO i=1,iim
	538	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
	539	ENDDO
	540	pdhfi(iip1,j,l) = pdhfi(1,j,l)
	541	ENDDO
	542
	543	ENDDO
	544
	545
	546	c 62. humidite specifique
	547	c ---------------------
[1279]	548	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	549	! DO iq=1,nqtot
	550	! DO l=1,llm
	551	! DO i=1,iip1
	552	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	553	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
	554	! ENDDO
	555	! DO j=2,jjm
	556	! ig0=1+(j-2)*iim
	557	! DO i=1,iim
	558	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
	559	! ENDDO
	560	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
	561	! ENDDO
	562	! ENDDO
	563	! ENDDO
[524]	564
	565	c 63. traceurs
	566	c ------------
	567	C initialisation des tendances
[1279]	568	pdqfi(:,:,:,:)=0.
[524]	569	C
[1146]	570	DO iq=1,nqtot
[524]	571	iiq=niadv(iq)
	572	DO l=1,llm
	573	DO i=1,iip1
	574	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
	575	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
	576	ENDDO
	577	DO j=2,jjm
	578	ig0=1+(j-2)*iim
	579	DO i=1,iim
	580	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
	581	ENDDO
	582	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
	583	ENDDO
	584	ENDDO
	585	ENDDO
	586
	587	c 65. champ u:
	588	c ------------
	589
	590	DO l=1,llm
	591
	592	DO i=1,iip1
	593	pdufi(i,1,l) = 0.
	594	pdufi(i,jjp1,l) = 0.
	595	ENDDO
	596
	597	DO j=2,jjm
	598	ig0=1+(j-2)*iim
	599	DO i=1,iim-1
	600	pdufi(i,j,l)=
	601	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
	602	ENDDO
	603	pdufi(iim,j,l)=
	604	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
	605	pdufi(iip1,j,l)=pdufi(1,j,l)
	606	ENDDO
	607
	608	ENDDO
	609
	610
	611	c 67. champ v:
	612	c ------------
	613
	614	DO l=1,llm
	615
	616	DO j=2,jjm-1
	617	ig0=1+(j-2)*iim
	618	DO i=1,iim
	619	pdvfi(i,j,l)=
	620	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
	621	ENDDO
	622	pdvfi(iip1,j,l) = pdvfi(1,j,l)
	623	ENDDO
	624	ENDDO
	625
	626
	627	c 68. champ v pres des poles:
	628	c ---------------------------
	629	c v = U * cos(long) + V * SIN(long)
	630
	631	DO l=1,llm
	632
	633	DO i=1,iim
	634	pdvfi(i,1,l)=
	635	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	636	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
	637	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
	638	pdvfi(i,1,l)=
	639	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	640	pdvfi(i,jjm,l)=
	641	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
	642	ENDDO
	643
	644	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	645	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	646
	647	ENDDO
	648
	649	c-----------------------------------------------------------------------
	650
	651	700 CONTINUE
	652
	653	firstcal = .FALSE.
	654
	655	RETURN
	656	END

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

Download in other formats:

Original Format