Context Navigation

calfis.F @ 1972

Last change on this file since 1972 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.4 KB

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

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

Download in other formats:

Original Format