Context Navigation

calfis.F @ 5434

Last change on this file since 5434 was 4482, checked in by lguez, 21 months ago
Sync latest trunk changes to branch LMDZ_ECRad
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: 21.3 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: calfis.F 4482 2023-03-29 13:14:27Z 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	$ pdufi,
	24	$ pdvfi,
	25	$ pdhfi,
	26	$ pdqfi,
	27	$ pdpsfi)
	28	c
	29	c Auteur : P. Le Van, F. Hourdin
	30	c .........
[4056]	31	USE infotrac, ONLY: nqtot, tracers
[1987]	32	USE control_mod, ONLY: planet_type, nsplit_phys
[2418]	33	#ifdef CPP_PHYS
	34	USE callphysiq_mod, ONLY: call_physiq
	35	#endif
[2597]	36	USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
[2600]	37	USE comvert_mod, ONLY: preff, presnivs
[2597]	38
[524]	39	IMPLICIT NONE
	40	c=======================================================================
	41	c
	42	c 1. rearrangement des tableaux et transformation
	43	c variables dynamiques > variables physiques
	44	c 2. calcul des termes physiques
	45	c 3. retransformation des tendances physiques en tendances dynamiques
	46	c
	47	c remarques:
	48	c ----------
	49	c
	50	c - les vents sont donnes dans la physique par leurs composantes
	51	c naturelles.
	52	c - la variable thermodynamique de la physique est une variable
	53	c intensive : T
	54	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
	55	c - les deux seules variables dependant de la geometrie necessaires
	56	c pour la physique sont la latitude pour le rayonnement et
	57	c l'aire de la maille quand on veut integrer une grandeur
	58	c horizontalement.
	59	c - les points de la physique sont les points scalaires de la
	60	c la dynamique; numerotation:
	61	c 1 pour le pole nord
	62	c (jjm-1)*iim pour l'interieur du domaine
	63	c ngridmx pour le pole sud
	64	c ---> ngridmx=2+(jjm-1)*iim
	65	c
	66	c Input :
	67	c -------
	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
[2600]	93	include "dimensions.h"
	94	include "paramet.h"
[524]	95
[1146]	96	INTEGER ngridmx
[524]	97	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	98
[2600]	99	include "comgeom2.h"
	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
[2418]	124	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
[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
	133
	134	c Local variables :
	135	c -----------------
	136
[4056]	137	INTEGER i,j,l,ig0,ig,iq,itr
[524]	138	REAL zpsrf(ngridmx)
	139	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
	140	REAL zphi(ngridmx,llm),zphis(ngridmx)
	141	c
[2333]	142	REAL zrot(iip1,jjm,llm) ! AdlC May 2014
[2604]	143	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
	144	REAL zrfi(ngridmx,llm) ! relative wind vorticity
[1146]	145	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
[2604]	146	REAL zpk(ngridmx,llm)
[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
[524]	173	c
	174	c-----------------------------------------------------------------------
	175	c
	176	c 1. Initialisations :
	177	c --------------------
	178	c
[1279]	179	c
	180	IF ( firstcal ) THEN
	181	debut = .TRUE.
	182	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
[1403]	183	write(lunout,*) 'STOP dans calfis'
	184	write(lunout,*)
	185	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
	186	write(lunout,*) ' ngridmx jjm iim '
	187	write(lunout,*) ngridmx,jjm,iim
[4482]	188	call abort_gcm("calfis", "", 1)
[1279]	189	ENDIF
[524]	190	ELSE
[1279]	191	debut = .FALSE.
	192	ENDIF ! of IF (firstcal)
[524]	193
	194	c
	195	c
	196	c-----------------------------------------------------------------------
	197	c 40. transformation des variables dynamiques en variables physiques:
	198	c ---------------------------------------------------------------
	199
	200	c 41. pressions au sol (en Pascals)
	201	c ----------------------------------
	202
	203
	204	zpsrf(1) = pps(1,1)
	205
	206	ig0 = 2
	207	DO j = 2,jjm
	208	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
	209	ig0 = ig0+iim
	210	ENDDO
	211
	212	zpsrf(ngridmx) = pps(1,jjp1)
	213
	214
[2604]	215	c 42. pression intercouches et fonction d'Exner:
[524]	216	c
	217	c -----------------------------------------------------------------
	218	c .... zplev definis aux (llm +1) interfaces des couches ....
	219	c .... zplay definis aux ( llm ) milieux des couches ....
	220	c -----------------------------------------------------------------
	221
	222	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
	223	c
	224	unskap = 1./ kappa
	225	c
[2604]	226	DO l = 1, llm
	227	zpk( 1,l ) = ppk(1,1,l)
[524]	228	zplev( 1,l ) = pp(1,1,l)
	229	ig0 = 2
	230	DO j = 2, jjm
	231	DO i =1, iim
[2604]	232	zpk( ig0,l ) = ppk(i,j,l)
[524]	233	zplev( ig0,l ) = pp(i,j,l)
	234	ig0 = ig0 +1
	235	ENDDO
	236	ENDDO
[2604]	237	zpk( ngridmx,l ) = ppk(1,jjp1,l)
[524]	238	zplev( ngridmx,l ) = pp(1,jjp1,l)
	239	ENDDO
[2604]	240	zplev( 1,llmp1 ) = pp(1,1,llmp1)
	241	ig0 = 2
	242	DO j = 2, jjm
	243	DO i =1, iim
	244	zplev( ig0,llmp1 ) = pp(i,j,llmp1)
	245	ig0 = ig0 +1
	246	ENDDO
	247	ENDDO
	248	zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
[524]	249	c
	250	c
	251
	252	c 43. temperature naturelle (en K) et pressions milieux couches .
	253	c ---------------------------------------------------------------
	254
	255	DO l=1,llm
	256
	257	pksurcp = ppk(1,1,l) / cpp
	258	zplay(1,l) = preff * pksurcp ** unskap
	259	ztfi(1,l) = pteta(1,1,l) * pksurcp
	260	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
	261	ig0 = 2
	262
	263	DO j = 2, jjm
	264	DO i = 1, iim
	265	pksurcp = ppk(i,j,l) / cpp
	266	zplay(ig0,l) = preff * pksurcp ** unskap
	267	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
	268	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
	269	ig0 = ig0 + 1
	270	ENDDO
	271	ENDDO
	272
	273	pksurcp = ppk(1,jjp1,l) / cpp
	274	zplay(ig0,l) = preff * pksurcp ** unskap
	275	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
	276	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
	277
	278	ENDDO
	279
	280	c 43.bis traceurs
	281	c ---------------
	282	c
[4056]	283	itr=0
[1146]	284	DO iq=1,nqtot
[4056]	285	IF(.NOT.tracers(iq)%isAdvected) CYCLE
	286	itr = itr + 1
[524]	287	DO l=1,llm
[4056]	288	zqfi(1,l,itr) = pq(1,1,l,iq)
	289	ig0 = 2
[524]	290	DO j=2,jjm
	291	DO i = 1, iim
[4056]	292	zqfi(ig0,l,itr) = pq(i,j,l,iq)
[524]	293	ig0 = ig0 + 1
	294	ENDDO
	295	ENDDO
[4056]	296	zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
[524]	297	ENDDO
	298	ENDDO
	299
	300	c convergence dynamique pour les traceurs "EAU"
[1279]	301	! Earth-specific treatment of first 2 tracers (water)
	302	if (planet_type=="earth") then
	303	DO iq=1,2
[524]	304	DO l=1,llm
	305	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
	306	ig0 = 2
	307	DO j=2,jjm
	308	DO i = 1, iim
	309	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
	310	ig0 = ig0 + 1
	311	ENDDO
	312	ENDDO
	313	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
	314	ENDDO
[1279]	315	ENDDO
	316	endif ! of if (planet_type=="earth")
[524]	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
[1403]	325	DO ig=1,ngridmx
	326	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	327	ENDDO
[524]	328	ENDDO
	329
	330	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
[960]	331	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
	332	c de masse est calclue dans advtrac.F
	333	c DO l=1,llm
	334	c pvervel(1,l)=pw(1,1,l) * g /apoln
	335	c ig0=2
	336	c DO j=2,jjm
	337	c DO i = 1, iim
	338	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
	339	c ig0 = ig0 + 1
	340	c ENDDO
	341	c ENDDO
	342	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
	343	c ENDDO
[524]	344
	345	c
	346	c 45. champ u:
	347	c ------------
	348
	349	DO 50 l=1,llm
	350
	351	DO 25 j=2,jjm
	352	ig0 = 1+(j-2)*iim
	353	zufi(ig0+1,l)= 0.5 *
	354	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
	355	pcvgu(ig0+1,l)= 0.5 *
	356	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
	357	DO 10 i=2,iim
	358	zufi(ig0+i,l)= 0.5 *
	359	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
	360	pcvgu(ig0+i,l)= 0.5 *
	361	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
	362	10 CONTINUE
	363	25 CONTINUE
	364
	365	50 CONTINUE
	366
	367
[2333]	368	C Alvaro de la Camara (May 2014)
	369	C 46.1 Calcul de la vorticite et passage sur la grille physique
	370	C --------------------------------------------------------------
	371	DO l=1,llm
	372	do i=1,iim
	373	do j=1,jjm
	374	zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
	375	$ + pucov(i,j+1,l) - pucov(i,j,l))
	376	$ / (cu(i,j)+cu(i,j+1))
	377	$ / (cv(i+1,j)+cv(i,j)) *4
	378	enddo
	379	enddo
	380	ENDDO
	381
[524]	382	c 46.champ v:
	383	c -----------
	384
	385	DO l=1,llm
	386	DO j=2,jjm
	387	ig0=1+(j-2)*iim
	388	DO i=1,iim
	389	zvfi(ig0+i,l)= 0.5 *
	390	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
	391	pcvgv(ig0+i,l)= 0.5 *
	392	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
	393	ENDDO
[2333]	394	zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
	395	& +zrot(1,j-1,l)+zrot(1,j,l))
	396	DO i=2,iim
	397	zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
	398	$ +zrot(i,j-1,l)+zrot(i,j,l)) ! AdlC MAY 2014
	399	ENDDO
[524]	400	ENDDO
	401	ENDDO
	402
	403
	404	c 47. champs de vents aux pole nord
	405	c ------------------------------
	406	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	407	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	408
	409	DO l=1,llm
	410
	411	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	412	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
	413	DO i=2,iim
	414	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	415	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
	416	ENDDO
	417
	418	DO i=1,iim
	419	zcos(i) = COS(rlonv(i))*z1(i)
	420	zcosbis(i)= COS(rlonv(i))*z1bis(i)
	421	zsin(i) = SIN(rlonv(i))*z1(i)
	422	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
	423	ENDDO
	424
	425	zufi(1,l) = SSUM(iim,zcos,1)/pi
	426	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
	427	zvfi(1,l) = SSUM(iim,zsin,1)/pi
	428	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
[2333]	429	zrfi(1, l) = 0.
[524]	430	ENDDO
	431
	432
	433	c 48. champs de vents aux pole sud:
	434	c ---------------------------------
	435	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	436	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	437
	438	DO l=1,llm
	439
	440	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	441	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
	442	DO i=2,iim
	443	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	444	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
[1403]	445	ENDDO
[524]	446
	447	DO i=1,iim
	448	zcos(i) = COS(rlonv(i))*z1(i)
	449	zcosbis(i) = COS(rlonv(i))*z1bis(i)
	450	zsin(i) = SIN(rlonv(i))*z1(i)
	451	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
[1403]	452	ENDDO
[524]	453
	454	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
	455	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
	456	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
	457	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
[2333]	458	zrfi(ngridmx, l) = 0.
[524]	459	ENDDO
[644]	460	c
[960]	461	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
[524]	462	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
	463
	464	c-----------------------------------------------------------------------
	465	c Appel de la physique:
	466	c ---------------------
	467
	468
[1403]	469
[1407]	470	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
[1403]	471	zdt_split=dtphys/nsplit_phys
	472	zdufic(:,:)=0.
	473	zdvfic(:,:)=0.
	474	zdtfic(:,:)=0.
	475	zdqfic(:,:,:)=0.
	476
[1615]	477	#ifdef CPP_PHYS
[1403]	478
[1615]	479	do isplit=1,nsplit_phys
	480
[1403]	481	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
	482	debut_split=debut.and.isplit==1
	483	lafin_split=lafin.and.isplit==nsplit_phys
	484
[4056]	485	! if (planet_type=="earth") then
[4046]	486	CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
[2418]	487	& debut_split,lafin_split,
	488	& jD_cur,jH_cur_split,zdt_split,
	489	& zplev,zplay,
[2604]	490	& zpk,zphi,zphis,
[2418]	491	& presnivs,
	492	& zufi,zvfi,zrfi,ztfi,zqfi,
	493	& flxwfi,pducov,
	494	& zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
	495	!
	496	! else if ( planet_type=="generic" ) then
	497	!
	498	! CALL physiq (ngridmx, !! ngrid
	499	! . llm, !! nlayer
	500	! . nqtot, !! nq
[4046]	501	! . tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
[2418]	502	! . debut_split, !! firstcall
	503	! . lafin_split, !! lastcall
	504	! . jD_cur, !! pday. see leapfrog
	505	! . jH_cur_split, !! ptime "fraction of day"
	506	! . zdt_split, !! ptimestep
	507	! . zplev, !! pplev
	508	! . zplay, !! pplay
	509	! . zphi, !! pphi
	510	! . zufi, !! pu
	511	! . zvfi, !! pv
	512	! . ztfi, !! pt
	513	! . zqfi, !! pq
	514	! . flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
	515	! . zdufi, !! pdu
	516	! . zdvfi, !! pdv
	517	! . zdtfi, !! pdt
	518	! . zdqfi, !! pdq
	519	! . zdpsrf, !! pdpsrf
	520	! . tracerdyn) !! tracerdyn <-- utilite ???
	521	!
	522	! endif ! of if (planet_type=="earth")
[1654]	523
[1403]	524	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
	525	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
	526	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
	527	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
	528
	529	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
	530	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
	531	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
	532	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
	533
[1615]	534	enddo ! of do isplit=1,nsplit_phys
	535
	536	#endif
	537	! of #ifdef CPP_PHYS
	538
[1403]	539	zdufi(:,:)=zdufic(:,:)/nsplit_phys
	540	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
	541	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
	542	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
	543
[524]	544
	545	500 CONTINUE
	546
	547	c-----------------------------------------------------------------------
	548	c transformation des tendances physiques en tendances dynamiques:
	549	c ---------------------------------------------------------------
	550
	551	c tendance sur la pression :
	552	c -----------------------------------
	553
	554	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
	555	c
	556	c 62. enthalpie potentielle
	557	c ---------------------
	558
	559	DO l=1,llm
	560
	561	DO i=1,iip1
	562	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	563	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
	564	ENDDO
	565
	566	DO j=2,jjm
	567	ig0=1+(j-2)*iim
	568	DO i=1,iim
	569	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
	570	ENDDO
	571	pdhfi(iip1,j,l) = pdhfi(1,j,l)
	572	ENDDO
	573
	574	ENDDO
	575
	576
	577	c 62. humidite specifique
	578	c ---------------------
[1279]	579	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	580	! DO iq=1,nqtot
	581	! DO l=1,llm
	582	! DO i=1,iip1
	583	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	584	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
	585	! ENDDO
	586	! DO j=2,jjm
	587	! ig0=1+(j-2)*iim
	588	! DO i=1,iim
	589	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
	590	! ENDDO
	591	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
	592	! ENDDO
	593	! ENDDO
	594	! ENDDO
[524]	595
	596	c 63. traceurs
	597	c ------------
	598	C initialisation des tendances
[1279]	599	pdqfi(:,:,:,:)=0.
[524]	600	C
[4056]	601	itr = 0
[1146]	602	DO iq=1,nqtot
[4056]	603	IF(.NOT.tracers(iq)%isAdvected) CYCLE
	604	itr = itr + 1
[524]	605	DO l=1,llm
	606	DO i=1,iip1
[4056]	607	pdqfi(i,1,l,iq) = zdqfi(1,l,itr)
	608	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
[524]	609	ENDDO
	610	DO j=2,jjm
	611	ig0=1+(j-2)*iim
	612	DO i=1,iim
[4056]	613	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
[524]	614	ENDDO
[4056]	615	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
[524]	616	ENDDO
	617	ENDDO
	618	ENDDO
	619
	620	c 65. champ u:
	621	c ------------
	622
	623	DO l=1,llm
	624
	625	DO i=1,iip1
	626	pdufi(i,1,l) = 0.
	627	pdufi(i,jjp1,l) = 0.
	628	ENDDO
	629
	630	DO j=2,jjm
	631	ig0=1+(j-2)*iim
	632	DO i=1,iim-1
	633	pdufi(i,j,l)=
	634	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
	635	ENDDO
	636	pdufi(iim,j,l)=
	637	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
	638	pdufi(iip1,j,l)=pdufi(1,j,l)
	639	ENDDO
	640
	641	ENDDO
	642
	643
	644	c 67. champ v:
	645	c ------------
	646
	647	DO l=1,llm
	648
	649	DO j=2,jjm-1
	650	ig0=1+(j-2)*iim
	651	DO i=1,iim
	652	pdvfi(i,j,l)=
	653	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
	654	ENDDO
	655	pdvfi(iip1,j,l) = pdvfi(1,j,l)
	656	ENDDO
	657	ENDDO
	658
	659
	660	c 68. champ v pres des poles:
	661	c ---------------------------
	662	c v = U * cos(long) + V * SIN(long)
	663
	664	DO l=1,llm
	665
	666	DO i=1,iim
	667	pdvfi(i,1,l)=
	668	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	669	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
	670	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
	671	pdvfi(i,1,l)=
	672	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	673	pdvfi(i,jjm,l)=
	674	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
	675	ENDDO
	676
	677	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	678	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	679
	680	ENDDO
	681
	682	c-----------------------------------------------------------------------
	683
	684	700 CONTINUE
	685
	686	firstcal = .FALSE.
	687
	688	RETURN
	689	END

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

Download in other formats:

Original Format