Context Navigation

calfis.F @ 3553

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

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

Download in other formats:

Original Format