Context Navigation

calfis.F @ 1384

Last change on this file since 1384 was 1312, checked in by emillour, 11 years ago
Generic and Mars GCM: LMDZ.COMMON dynamics sends mass flux to physics and not vertical velocity. Add the computation of vertical velocity from input mass flux in the physics, and also modify "old" LMDZ.GENERIC and LMDZ.MARS dynamics to be consistent. EM
File size: 29.5 KB

Rev	Line
[1]	1	!
	2	! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
	3	!
	4	C
	5	C
	6	SUBROUTINE calfis(lafin,
	7	$ jD_cur, jH_cur,
	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 .........
[1189]	31	USE infotrac, ONLY: nqtot, niadv, tname
	32	USE control_mod, ONLY: planet_type, nsplit_phys
[1056]	33	USE write_field
[1017]	34	USE cpdet_mod, only: t2tpot,tpot2t
[1056]	35
	36	! used only for zonal averages
	37	USE moyzon_mod
[1]	38
	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
[5]	80	c pdhfi tendency for the potential temperature (K/s)
[1]	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"
[1056]	96	#include "logic.h"
[1]	97
	98	INTEGER ngridmx
	99	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	100
	101	#include "comconst.h"
	102	#include "comvert.h"
	103	#include "comgeom2.h"
	104	#include "iniprint.h"
	105
	106	c Arguments :
	107	c -----------
[1189]	108	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
	109	REAL,INTENT(IN) :: jD_cur, jH_cur
	110	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
	111	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
	112	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
	113	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
	114	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
	115	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
	116	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
[1]	117
[1189]	118	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
	119	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
	120	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
[6]	121	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
	122	! qui lui meme ne sert a rien dans la routine telle qu'elle est
	123	! ecrite, et que j'ai donc commente....
[1189]	124	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
	125	! NB: pdq is only used to compute pcvgq which is in fact not used...
[1056]	126
[1189]	127	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
	128	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
	129	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
[1312]	130	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
[1056]	131
[1189]	132	! tendencies (in */s) from the physics
	133	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
	134	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
	135	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
	136	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
	137	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
[1]	138
	139
	140	c Local variables :
	141	c -----------------
	142
	143	INTEGER i,j,l,ig0,ig,iq,iiq
	144	REAL zpsrf(ngridmx)
	145	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
	146	REAL zphi(ngridmx,llm),zphis(ngridmx)
[1056]	147
[1]	148	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
	149	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
[5]	150	! ADAPTATION GCM POUR CP(T)
	151	REAL zteta(ngridmx,llm)
	152	REAL zpk(ngridmx,llm)
[1056]	153
[5]	154	! RQ SL 13/10/10:
	155	! Ces calculs ne servent pas.
	156	! Si necessaire, decommenter ces variables et les calculs...
	157	! REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
	158	! REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
[1056]	159
[1]	160	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
	161	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
	162	REAL zdpsrf(ngridmx)
[1056]	163
[1]	164	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
	165	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
	166	REAL jH_cur_split,zdt_split
	167	LOGICAL debut_split,lafin_split
	168	INTEGER isplit
	169
	170	REAL zsin(iim),zcos(iim),z1(iim)
	171	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	172	REAL unskap, pksurcp
[1056]	173	save unskap
	174
[1]	175	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
	176
	177	REAL SSUM
	178
[1189]	179	LOGICAL,SAVE :: firstcal=.true., debut=.true.
[1]	180	! REAL rdayvrai
[841]	181
	182	LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
[1056]	183
	184	! For Titan only right now:
	185	! to allow for 2D computation of microphys and chemistry
	186	LOGICAL,save :: flag_moyzon
[1076]	187	REAL,allocatable,save :: tmpvar(:,:)
	188	REAL,allocatable,save :: tmpvarp1(:,:)
	189	REAL,allocatable,save :: tmpvarbar(:)
	190	REAL,allocatable,save :: tmpvarbarp1(:)
[1126]	191	real :: zz1,zz2
[1056]	192
[1]	193	c-----------------------------------------------------------------------
[1056]	194
[1]	195	c 1. Initialisations :
	196	c --------------------
[1056]	197
	198
[1]	199	IF ( firstcal ) THEN
	200	debut = .TRUE.
	201	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
	202	write(lunout,*) 'STOP dans calfis'
	203	write(lunout,*)
	204	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
	205	write(lunout,*) ' ngridmx jjm iim '
	206	write(lunout,*) ngridmx,jjm,iim
	207	STOP
	208	ENDIF
[1056]	209
	210	unskap = 1./ kappa
	211
	212	flag_moyzon = .false.
	213	if(moyzon_ch.or.moyzon_mu) then
	214	flag_moyzon = .true.
[1071]	215	allocate(tmpvar(iip1,llm))
	216	allocate(tmpvarp1(iip1,llmp1))
	217	allocate(tmpvarbar(llm))
	218	allocate(tmpvarbarp1(llmp1))
[1056]	219	endif
[1071]	220
[1256]	221	if (flag_moyzon) call moyzon_init(ngridmx,llm,nqtot)
[1056]	222
	223	c----------------------------------------------
	224	c moyennes globales pour le profil de pression
[1126]	225	if(planet_type.eq."titan".or.planet_type.eq."venus") then
[1056]	226	ALLOCATE(plevmoy(llm+1))
	227	ALLOCATE(playmoy(llm))
	228	ALLOCATE(tmoy(llm))
	229	ALLOCATE(tetamoy(llm))
	230	ALLOCATE(pkmoy(llm))
[1126]	231	ALLOCATE(phimoy(0:llm))
	232	ALLOCATE(zlevmoy(llm+1))
	233	ALLOCATE(zlaymoy(llm))
[1056]	234	plevmoy=0.
	235	do l=1,llmp1
	236	do i=1,iip1
	237	do j=1,jjp1
	238	plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
	239	enddo
	240	enddo
	241	enddo
	242	tetamoy=0.
	243	pkmoy=0.
[1126]	244	phimoy=0.
	245	do i=1,iip1
	246	do j=1,jjp1
	247	phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
	248	enddo
	249	enddo
[1056]	250	do l=1,llm
	251	do i=1,iip1
	252	do j=1,jjp1
	253	tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
	254	pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
[1126]	255	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
[1056]	256	enddo
	257	enddo
	258	enddo
[1126]	259	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
[1056]	260	call tpot2t(llm,tetamoy,tmoy,pkmoy)
[1126]	261	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
	262	zlaymoy(:) = gradrad/(g*rad-phimoy(:))-rad
	263	zlevmoy(1) = phimoy(0)/g
	264	DO l=2,llm
	265	zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
	266	zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l))
	267	zlevmoy(l)=(zz1zlaymoy(l-1)+zz2zlaymoy(l))/(zz1+zz2)
	268	ENDDO
	269	zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
[1071]	270	c-------------------
[1056]	271	c + lat index
[1071]	272	allocate(klat(ngridmx))
	273	klat=0
	274	klat(1) = 1
	275	ig0 = 2
	276	DO j = 2,jjm
[1126]	277	do i=0,iim-1
	278	klat(ig0+i) = j
	279	enddo
	280	ig0 = ig0+iim
[1071]	281	ENDDO
	282	klat(ngridmx) = jjp1
[1126]	283	endif ! planet_type=titan
[1056]	284	c----------------------------------------------
[1]	285	ELSE
	286	debut = .FALSE.
	287	ENDIF ! of IF (firstcal)
	288
[1056]	289
[1]	290	c-----------------------------------------------------------------------
	291	c 40. transformation des variables dynamiques en variables physiques:
	292	c ---------------------------------------------------------------
	293
	294	c 41. pressions au sol (en Pascals)
	295	c ----------------------------------
	296
	297	zpsrf(1) = pps(1,1)
	298
	299	ig0 = 2
	300	DO j = 2,jjm
	301	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
	302	ig0 = ig0+iim
	303	ENDDO
	304
	305	zpsrf(ngridmx) = pps(1,jjp1)
	306
[1056]	307	c 42. pression intercouches et fonction d'Exner:
[1]	308
	309	c -----------------------------------------------------------------
	310	c .... zplev definis aux (llm +1) interfaces des couches ....
	311	c .... zplay definis aux ( llm ) milieux des couches ....
	312	c -----------------------------------------------------------------
	313
	314	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
[1056]	315
[5]	316	! ADAPTATION GCM POUR CP(T)
	317	DO l = 1, llm
	318	zpk( 1,l ) = ppk(1,1,l)
	319	zteta( 1,l ) = pteta(1,1,l)
[1]	320	zplev( 1,l ) = pp(1,1,l)
	321	ig0 = 2
	322	DO j = 2, jjm
	323	DO i =1, iim
[5]	324	zpk( ig0,l ) = ppk(i,j,l)
	325	zteta( ig0,l ) = pteta(i,j,l)
[1]	326	zplev( ig0,l ) = pp(i,j,l)
	327	ig0 = ig0 +1
	328	ENDDO
	329	ENDDO
[5]	330	zpk( ngridmx,l ) = ppk(1,jjp1,l)
	331	zteta( ngridmx,l ) = pteta(1,jjp1,l)
[1]	332	zplev( ngridmx,l ) = pp(1,jjp1,l)
	333	ENDDO
[5]	334	zplev( 1,llmp1 ) = pp(1,1,llmp1)
	335	ig0 = 2
	336	DO j = 2, jjm
	337	DO i =1, iim
	338	zplev( ig0,llmp1 ) = pp(i,j,llmp1)
	339	ig0 = ig0 +1
	340	ENDDO
	341	ENDDO
	342	zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
[1]	343
[1056]	344	if (flag_moyzon) then
	345	tmpvarp1(:,:) = pp(:,1,:)
	346	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
	347	zplevbar(1,:) = tmpvarbarp1
	348	tmpvar(:,:) = ppk(:,1,:)
	349	call moyzon(llm,tmpvar,tmpvarbar)
	350	zpkbar(1,:) = tmpvarbar
	351	tmpvar(:,:) = pteta(:,1,:)
	352	call moyzon(llm,tmpvar,tmpvarbar)
	353	ztetabar(1,:) = tmpvarbar
	354	call tpot2t(llm,ztetabar(1,:),ztfibar(1,:),zpkbar(1,:))
	355	ig0 = 2
	356	do j = 2, jjm
	357	tmpvarp1(:,:) = pp(:,j,:)
	358	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
	359	zplevbar(ig0,:) = tmpvarbarp1
	360	tmpvar(:,:) = ppk(:,j,:)
	361	call moyzon(llm,tmpvar,tmpvarbar)
	362	zpkbar(ig0,:) = tmpvarbar
	363	tmpvar(:,:) = pteta(:,j,:)
	364	call moyzon(llm,tmpvar,tmpvarbar)
	365	ztetabar(ig0,:) = tmpvarbar
	366	call tpot2t(llm,ztetabar(ig0,:),ztfibar(ig0,:),zpkbar(ig0,:))
	367	ig0 = ig0+1
	368	do i=2,iim
	369	zplevbar(ig0,:) = zplevbar(ig0-1,:)
	370	zpkbar(ig0,:) = zpkbar(ig0-1,:)
	371	ztetabar(ig0,:) = ztetabar(ig0-1,:)
	372	ztfibar(ig0,:) = ztfibar(ig0-1,:)
	373	ig0 = ig0+1
	374	enddo
	375	enddo
	376	tmpvarp1(:,:) = pp(:,jjp1,:)
	377	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
	378	zplevbar(ngridmx,:) = tmpvarbarp1
	379	tmpvar(:,:) = ppk(:,jjp1,:)
	380	call moyzon(llm,tmpvar,tmpvarbar)
	381	zpkbar(ngridmx,:) = tmpvarbar
	382	tmpvar(:,:) = pteta(:,jjp1,:)
	383	call moyzon(llm,tmpvar,tmpvarbar)
	384	ztetabar(ngridmx,:) = tmpvarbar
	385	call tpot2t(llm,ztetabar(ngridmx,:),
	386	. ztfibar(ngridmx,:),zpkbar(ngridmx,:))
	387	endif
	388
[1]	389	c 43. temperature naturelle (en K) et pressions milieux couches .
	390	c ---------------------------------------------------------------
	391
[5]	392	! ADAPTATION GCM POUR CP(T)
	393	call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
	394
[1]	395	DO l=1,llm
	396
	397	pksurcp = ppk(1,1,l) / cpp
	398	zplay(1,l) = preff * pksurcp ** unskap
[5]	399	! pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
[1]	400	ig0 = 2
	401
	402	DO j = 2, jjm
	403	DO i = 1, iim
	404	pksurcp = ppk(i,j,l) / cpp
	405	zplay(ig0,l) = preff * pksurcp ** unskap
[5]	406	! pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
[1]	407	ig0 = ig0 + 1
	408	ENDDO
	409	ENDDO
	410
	411	pksurcp = ppk(1,jjp1,l) / cpp
	412	zplay(ig0,l) = preff * pksurcp ** unskap
[5]	413	! pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
[1]	414
	415	ENDDO
	416
[1056]	417	if (flag_moyzon) then
	418	zplaybar(:,:) = preff * (zpkbar(:,:)/cpp)**unskap
	419	endif
	420
[5]	421	c 43.bis traceurs (tous intensifs)
[1]	422	c ---------------
[1056]	423
[1]	424	DO iq=1,nqtot
	425	DO l=1,llm
[1056]	426	zqfi(1,l,iq) = pq(1,1,l,iq)
[1]	427	ig0 = 2
	428	DO j=2,jjm
	429	DO i = 1, iim
[1056]	430	zqfi(ig0,l,iq) = pq(i,j,l,iq)
[1]	431	ig0 = ig0 + 1
	432	ENDDO
	433	ENDDO
[1056]	434	zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
[1]	435	ENDDO
[5]	436	ENDDO ! boucle sur traceurs
[1]	437
[1056]	438	if (flag_moyzon) then
	439	DO iq=1,nqtot
	440	! RQ: REVOIR A QUOI CA SERT... ET VERIFIER...
	441	! iiq=niadv(iq)
	442	! en fait, iiq=iq...
	443	! FIN RQ
	444	tmpvar(:,:) = pq(:,1,:,iq)
	445	call moyzon(llm,tmpvar,tmpvarbar)
	446	zqfibar(1,:,iq) = tmpvarbar
	447	ig0 = 2
	448	do j = 2, jjm
	449	tmpvar(:,:) = pq(:,j,:,iq)
	450	call moyzon(llm,tmpvar,tmpvarbar)
	451	zqfibar(ig0,:,iq) = tmpvarbar
	452	ig0 = ig0+1
	453	do i=2,iim
	454	zqfibar(ig0,:,iq) = zqfibar(ig0-1,:,iq)
	455	ig0 = ig0+1
	456	enddo
	457	enddo
	458	tmpvar(:,:) = pq(:,jjp1,:,iq)
	459	call moyzon(llm,tmpvar,tmpvarbar)
	460	zqfibar(ngridmx,:,iq) = tmpvarbar
	461	ENDDO ! of DO iq=1,nqtot
	462	endif
	463
	464	! DEBUG
	465	! do ig0=1,ngridmx
	466	! write(*,'(6(e13.5,1x))') zqfibar(ig0,1,10),zqfi(ig0,1,10),
	467	! . zqfibar(ig0,llm/2,10),zqfi(ig0,llm/2,10),
	468	! . zqfibar(ig0,llm,10),zqfi(ig0,llm,10)
	469	! enddo
	470	! stop
	471
[5]	472	!-----------------
	473	! RQ SL 13/10/10:
	474	! Ces calculs ne servent pas.
	475	! Si necessaire, decommenter ces variables et les calculs...
	476	!
	477	! convergence dynamique pour les traceurs "EAU"
[1]	478	! Earth-specific treatment of first 2 tracers (water)
[5]	479	! if (planet_type=="earth") then
	480	! DO iq=1,2
	481	! DO l=1,llm
	482	! pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
	483	! ig0 = 2
	484	! DO j=2,jjm
	485	! DO i = 1, iim
	486	! pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
	487	! ig0 = ig0 + 1
	488	! ENDDO
	489	! ENDDO
	490	! pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
	491	! ENDDO
	492	! ENDDO
	493	! endif ! of if (planet_type=="earth")
	494	!----------------
[1]	495
	496	c Geopotentiel calcule par rapport a la surface locale:
	497	c -----------------------------------------------------
	498
	499	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
	500	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
	501	DO l=1,llm
	502	DO ig=1,ngridmx
	503	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	504	ENDDO
	505	ENDDO
	506
[1056]	507	if (flag_moyzon) then
	508	tmpvar(:,1) = pphis(:,1)
	509	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
	510	zphisbar(1) = tmpvarbar(1)
	511	tmpvar(:,:) = pphi(:,1,:)
	512	call moyzon(llm,tmpvar,tmpvarbar)
	513	zphibar(1,:) = tmpvarbar
	514	ig0 = 2
	515	do j = 2, jjm
	516	tmpvar(:,1) = pphis(:,j)
	517	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
	518	zphisbar(ig0) = tmpvarbar(1)
	519	tmpvar(:,:) = pphi(:,j,:)
	520	call moyzon(llm,tmpvar,tmpvarbar)
	521	zphibar(ig0,:) = tmpvarbar
	522	ig0 = ig0+1
	523	do i=2,iim
	524	zphisbar(ig0) = zphisbar(ig0-1)
	525	zphibar(ig0,:) = zphibar(ig0-1,:)
	526	ig0 = ig0+1
	527	enddo
	528	enddo
	529	tmpvar(:,1) = pphis(:,jjp1)
	530	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
	531	zphisbar(ngridmx) = tmpvarbar(1)
	532	tmpvar(:,:) = pphi(:,jjp1,:)
	533	call moyzon(llm,tmpvar,tmpvarbar)
	534	zphibar(ngridmx,:) = tmpvarbar
	535	endif
	536
[1]	537	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
	538	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
	539	c de masse est calclue dans advtrac.F
	540	c DO l=1,llm
	541	c pvervel(1,l)=pw(1,1,l) * g /apoln
	542	c ig0=2
	543	c DO j=2,jjm
	544	c DO i = 1, iim
	545	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
	546	c ig0 = ig0 + 1
	547	c ENDDO
	548	c ENDDO
	549	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
	550	c ENDDO
	551
	552	c
	553	c 45. champ u:
	554	c ------------
	555
	556	DO 50 l=1,llm
	557
	558	DO 25 j=2,jjm
	559	ig0 = 1+(j-2)*iim
	560	zufi(ig0+1,l)= 0.5 *
	561	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
[5]	562	! pcvgu(ig0+1,l)= 0.5 *
	563	! $ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
[1]	564	DO 10 i=2,iim
	565	zufi(ig0+i,l)= 0.5 *
	566	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
[5]	567	! pcvgu(ig0+i,l)= 0.5 *
	568	! $ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
[1]	569	10 CONTINUE
	570	25 CONTINUE
	571
	572	50 CONTINUE
	573
	574
	575	c 46.champ v:
	576	c -----------
	577
	578	DO l=1,llm
	579	DO j=2,jjm
	580	ig0=1+(j-2)*iim
	581	DO i=1,iim
	582	zvfi(ig0+i,l)= 0.5 *
	583	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
[5]	584	c pcvgv(ig0+i,l)= 0.5 *
	585	c $ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
[1]	586	ENDDO
	587	ENDDO
	588	ENDDO
	589
	590
	591	c 47. champs de vents aux pole nord
	592	c ------------------------------
	593	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	594	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	595
	596	DO l=1,llm
	597
	598	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	599	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
	600	DO i=2,iim
	601	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	602	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
	603	ENDDO
	604
	605	DO i=1,iim
	606	zcos(i) = COS(rlonv(i))*z1(i)
	607	zcosbis(i)= COS(rlonv(i))*z1bis(i)
	608	zsin(i) = SIN(rlonv(i))*z1(i)
	609	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
	610	ENDDO
	611
	612	zufi(1,l) = SSUM(iim,zcos,1)/pi
[5]	613	! pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
[1]	614	zvfi(1,l) = SSUM(iim,zsin,1)/pi
[5]	615	! pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
[1]	616
	617	ENDDO
	618
	619
	620	c 48. champs de vents aux pole sud:
	621	c ---------------------------------
	622	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	623	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	624
	625	DO l=1,llm
	626
	627	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	628	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
	629	DO i=2,iim
	630	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	631	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
	632	ENDDO
	633
	634	DO i=1,iim
	635	zcos(i) = COS(rlonv(i))*z1(i)
	636	zcosbis(i) = COS(rlonv(i))*z1bis(i)
	637	zsin(i) = SIN(rlonv(i))*z1(i)
	638	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
	639	ENDDO
	640
	641	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
[5]	642	! pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
[1]	643	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
[5]	644	! pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
[1]	645
	646	ENDDO
	647	c
	648	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
	649	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
	650
	651	c-----------------------------------------------------------------------
	652	c Appel de la physique:
	653	c ---------------------
	654
[6]	655	! Appel de la physique: pose probleme quand on tourne
	656	! SANS physique, car physiq.F est dans le repertoire phy[]...
	657	! Il faut une cle CPP_PHYS
[1]	658
[6]	659	! Le fait que les arguments de physiq soient differents selon les planetes
	660	! ne pose pas de probleme a priori.
	661
[1]	662	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
	663	zdt_split=dtphys/nsplit_phys
	664	zdufic(:,:)=0.
	665	zdvfic(:,:)=0.
	666	zdtfic(:,:)=0.
	667	zdqfic(:,:,:)=0.
	668
[776]	669	#ifdef CPP_PHYS
	670
[1]	671	do isplit=1,nsplit_phys
	672
	673	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
	674	debut_split=debut.and.isplit==1
	675	lafin_split=lafin.and.isplit==nsplit_phys
	676
[6]	677	if (planet_type.eq."earth") then
[1]	678	CALL physiq (ngridmx,
	679	. llm,
	680	. debut_split,
	681	. lafin_split,
	682	. jD_cur,
	683	. jH_cur_split,
	684	. zdt_split,
	685	. zplev,
	686	. zplay,
	687	. zphi,
	688	. zphis,
	689	. presnivs,
	690	. zufi,
	691	. zvfi,
	692	. ztfi,
	693	. zqfi,
	694	. flxwfi,
	695	. zdufi,
	696	. zdvfi,
	697	. zdtfi,
	698	. zdqfi,
	699	. zdpsrf,
[1302]	700	. pducov)
[1]	701
[841]	702	else if ( planet_type=="generic" ) then
	703
	704	CALL physiq (ngridmx, !! ngrid
	705	. llm, !! nlayer
	706	. nqtot, !! nq
	707	. tname, !! tracer names from dynamical core (given in infotrac)
	708	. debut_split, !! firstcall
	709	. lafin_split, !! lastcall
	710	. jD_cur, !! pday. see leapfrog
	711	. jH_cur_split, !! ptime "fraction of day"
	712	. zdt_split, !! ptimestep
	713	. zplev, !! pplev
	714	. zplay, !! pplay
	715	. zphi, !! pphi
	716	. zufi, !! pu
	717	. zvfi, !! pv
	718	. ztfi, !! pt
	719	. zqfi, !! pq
	720	. flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
	721	. zdufi, !! pdu
	722	. zdvfi, !! pdv
	723	. zdtfi, !! pdt
	724	. zdqfi, !! pdq
	725	. zdpsrf, !! pdpsrf
	726	. tracerdyn) !! tracerdyn <-- utilite ???
	727
[1107]	728	else if ( planet_type=="mars" ) then
[1]	729
[1107]	730	CALL physiq (ngridmx, ! ngrid
	731	. llm, ! nlayer
	732	. nqtot, ! nq
	733	. debut_split, ! firstcall
	734	. lafin_split, ! lastcall
	735	. jD_cur, ! pday
	736	. jH_cur_split, ! ptime
	737	. zdt_split, ! ptimestep
	738	. zplev, ! pplev
	739	. zplay, ! pplay
	740	. zphi, ! pphi
	741	. zufi, ! pu
	742	. zvfi, ! pv
	743	. ztfi, ! pt
	744	. zqfi, ! pq
	745	. flxwfi, ! pw
	746	. zdufi, ! pdu
	747	. zdvfi, ! pdv
	748	. zdtfi, ! pdt
	749	. zdqfi, ! pdq
	750	. zdpsrf, ! pdpsrf
	751	. tracerdyn) ! tracerdyn (somewhat obsolete)
[5]	752
[1129]	753	else if ((planet_type=="titan").or.(planet_type=="venus")) then
	754
	755	CALL physiq (ngridmx,
	756	. llm,
	757	. nqtot,
	758	. debut_split,
	759	. lafin_split,
	760	. jD_cur,
	761	. jH_cur_split,
	762	. zdt_split,
	763	. zplev,
	764	. zplay,
	765	. zpk,
	766	. zphi,
	767	. zphis,
	768	. presnivs,
	769	. zufi,
	770	. zvfi,
	771	. ztfi,
	772	. zqfi,
	773	. flxwfi,
	774	. zdufi,
	775	. zdvfi,
	776	. zdtfi,
	777	. zdqfi,
	778	. zdpsrf)
	779
[1107]	780	else ! unknown "planet_type"
	781
	782	write(lunout,*) "calfis_p: error, unknown planet_type: ",
	783	& trim(planet_type)
	784	stop
	785
[6]	786	endif ! planet_type
[1]	787
[6]	788	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
	789	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
	790	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
	791	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
	792
	793	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
	794	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
	795	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
	796	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
	797
[776]	798	enddo ! of do isplit=1,nsplit_phys
	799
[1056]	800	! ATTENTION...
	801	if (flag_moyzon.and.(nsplit_phys.ne.1)) then
	802	print*,"WARNING ! flag_moyzon + nsplit_phys"
	803	print*,"zqfibar n'est pas implemente au cours des iterations"
	804	print*,"Donc a revoir..."
	805	stop
	806	endif
	807
[776]	808	#endif
	809	! #endif of #ifdef CPP_PHYS
	810
[6]	811	zdufi(:,:)=zdufic(:,:)/nsplit_phys
	812	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
	813	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
	814	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
	815
	816
[1]	817	500 CONTINUE
	818
	819	c-----------------------------------------------------------------------
	820	c transformation des tendances physiques en tendances dynamiques:
	821	c ---------------------------------------------------------------
	822
	823	c tendance sur la pression :
	824	c -----------------------------------
	825
	826	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
	827	c
	828	c 62. enthalpie potentielle
	829	c ---------------------
	830
[5]	831	! ADAPTATION GCM POUR CP(T)
	832	call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
[1]	833
	834	DO i=1,iip1
[5]	835	pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
	836	pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
[1]	837	ENDDO
	838
	839	DO j=2,jjm
	840	ig0=1+(j-2)*iim
	841	DO i=1,iim
[5]	842	pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
[1]	843	ENDDO
[5]	844	pdhfi(iip1,j,:) = pdhfi(1,j,:)
[1]	845	ENDDO
	846
	847
	848	c 62. humidite specifique
	849	c ---------------------
	850	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	851	! DO iq=1,nqtot
	852	! DO l=1,llm
	853	! DO i=1,iip1
	854	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	855	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
	856	! ENDDO
	857	! DO j=2,jjm
	858	! ig0=1+(j-2)*iim
	859	! DO i=1,iim
	860	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
	861	! ENDDO
	862	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
	863	! ENDDO
	864	! ENDDO
	865	! ENDDO
	866
[5]	867	c 63. traceurs (tous en intensifs)
[1]	868	c ------------
	869	C initialisation des tendances
	870	pdqfi(:,:,:,:)=0.
	871	C
[6]	872	DO iq=1,nqtot
[1]	873	iiq=niadv(iq)
	874	DO l=1,llm
	875	DO i=1,iip1
	876	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
	877	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
	878	ENDDO
	879	DO j=2,jjm
	880	ig0=1+(j-2)*iim
	881	DO i=1,iim
	882	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
	883	ENDDO
	884	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
	885	ENDDO
	886	ENDDO
[6]	887	ENDDO
[1]	888
	889	c 65. champ u:
	890	c ------------
	891
	892	DO l=1,llm
	893
	894	DO i=1,iip1
	895	pdufi(i,1,l) = 0.
	896	pdufi(i,jjp1,l) = 0.
	897	ENDDO
	898
	899	DO j=2,jjm
	900	ig0=1+(j-2)*iim
	901	DO i=1,iim-1
	902	pdufi(i,j,l)=
	903	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
	904	ENDDO
	905	pdufi(iim,j,l)=
	906	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
	907	pdufi(iip1,j,l)=pdufi(1,j,l)
	908	ENDDO
	909
	910	ENDDO
	911
	912
	913	c 67. champ v:
	914	c ------------
	915
	916	DO l=1,llm
	917
	918	DO j=2,jjm-1
	919	ig0=1+(j-2)*iim
	920	DO i=1,iim
	921	pdvfi(i,j,l)=
	922	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
	923	ENDDO
	924	pdvfi(iip1,j,l) = pdvfi(1,j,l)
	925	ENDDO
	926	ENDDO
	927
	928
	929	c 68. champ v pres des poles:
	930	c ---------------------------
	931	c v = U * cos(long) + V * SIN(long)
	932
	933	DO l=1,llm
	934
	935	DO i=1,iim
	936	pdvfi(i,1,l)=
	937	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	938	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
	939	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
	940	pdvfi(i,1,l)=
	941	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	942	pdvfi(i,jjm,l)=
	943	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
	944	ENDDO
	945
	946	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	947	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	948
	949	ENDDO
	950
	951	c-----------------------------------------------------------------------
	952
	953	700 CONTINUE
	954
	955	firstcal = .FALSE.
	956
	957	RETURN
	958	END

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

Download in other formats:

Original Format