Context Navigation

calfis_p.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: 42.4 KB

Rev	Line
[1]	1	!
	2	! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fairhead $
	3	!
	4	C
	5	C
	6	SUBROUTINE calfis_p(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)
[8]	28	#ifdef CPP_PHYS
[776]	29	! Ehouarn: if using (parallelized) physics
[1]	30	USE dimphy
[1543]	31	USE mod_phys_lmdz_mpi_data, mpi_root_xx=>mpi_master
	32	USE mod_phys_lmdz_omp_data, ONLY: klon_omp, klon_omp_begin
	33	USE mod_const_mpi, ONLY: COMM_LMDZ
[1]	34	USE mod_interface_dyn_phys
[1107]	35	! USE IOPHY
[1]	36	#endif
[1403]	37	#ifdef CPP_PARA
[1549]	38	USE parallel_lmdz, ONLY: omp_chunk, using_mpi, AllGather_Field,
	39	& jjb_u,jje_u,jjb_v,jje_v,
	40	& jj_begin_dyn=>jj_begin,jj_end_dyn=>jj_end
[1]	41	USE Write_Field
	42	Use Write_field_p
	43	USE Times
[1086]	44	USE cpdet_mod, only: tpot2t_p, t2tpot_p
[1056]	45	! used only for zonal averages
	46	USE moyzon_mod
[1403]	47	#endif
	48	USE infotrac, ONLY: nqtot, niadv, tname
	49	USE control_mod, ONLY: planet_type, nsplit_phys
[1422]	50	USE comvert_mod, ONLY: preff,presnivs
	51	USE comconst_mod, ONLY: daysec,dtvr,dtphys,kappa,cpp,g,rad,pi
	52	USE logic_mod, ONLY: moyzon_ch,moyzon_mu
[1549]	53	#ifdef CPP_PHYS
	54	USE callphysiq_mod, ONLY: call_physiq
	55	#endif
[1056]	56
[1]	57	IMPLICIT NONE
	58	c=======================================================================
	59	c
	60	c 1. rearrangement des tableaux et transformation
	61	c variables dynamiques > variables physiques
	62	c 2. calcul des termes physiques
	63	c 3. retransformation des tendances physiques en tendances dynamiques
	64	c
	65	c remarques:
	66	c ----------
	67	c
	68	c - les vents sont donnes dans la physique par leurs composantes
	69	c naturelles.
	70	c - la variable thermodynamique de la physique est une variable
	71	c intensive : T
	72	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
	73	c - les deux seules variables dependant de la geometrie necessaires
	74	c pour la physique sont la latitude pour le rayonnement et
	75	c l'aire de la maille quand on veut integrer une grandeur
	76	c horizontalement.
	77	c - les points de la physique sont les points scalaires de la
	78	c la dynamique; numerotation:
	79	c 1 pour le pole nord
	80	c (jjm-1)*iim pour l'interieur du domaine
	81	c ngridmx pour le pole sud
	82	c ---> ngridmx=2+(jjm-1)*iim
	83	c
	84	c Input :
	85	c -------
	86	c ecritphy frequence d'ecriture (en jours)de histphy
	87	c pucov covariant zonal velocity
	88	c pvcov covariant meridional velocity
	89	c pteta potential temperature
	90	c pps surface pressure
	91	c pmasse masse d'air dans chaque maille
	92	c pts surface temperature (K)
	93	c callrad clef d'appel au rayonnement
	94	c
	95	c Output :
	96	c --------
	97	c pdufi tendency for the natural zonal velocity (ms-1)
	98	c pdvfi tendency for the natural meridional velocity
[8]	99	c pdhfi tendency for the potential temperature (K/s)
[1]	100	c pdtsfi tendency for the surface temperature
	101	c
	102	c pdtrad radiative tendencies \ both input
	103	c pfluxrad radiative fluxes / and output
	104	c
	105	c=======================================================================
	106	c
	107	c-----------------------------------------------------------------------
	108	c
	109	c 0. Declarations :
	110	c ------------------
	111
[1459]	112	include "dimensions.h"
	113	include "paramet.h"
[1]	114
	115	INTEGER ngridmx
	116	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	117
[1459]	118	include "comgeom2.h"
	119	include "iniprint.h"
[1]	120	#ifdef CPP_MPI
	121	include 'mpif.h'
	122	#endif
	123	c Arguments :
	124	c -----------
[1189]	125	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
	126	REAL,INTENT(IN) :: jD_cur, jH_cur
	127	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
	128	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
	129	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
	130	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
	131	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
	132	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
	133	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
[1056]	134
[1189]	135	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
	136	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
	137	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
[8]	138	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
	139	! qui lui meme ne sert a rien dans la routine telle qu'elle est
	140	! ecrite, et que j'ai donc commente....
[1189]	141	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
	142	! NB: pdq is only used to compute pcvgq which is in fact not used...
[1056]	143
[1189]	144	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
	145	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
	146	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
[1312]	147	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
[1056]	148
[1189]	149	! tendencies (in */s) from the physics
	150	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
	151	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
	152	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
	153	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
	154	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
[1]	155
[1403]	156	#ifdef CPP_PARA
[8]	157	#ifdef CPP_PHYS
[1]	158	c Local variables :
	159	c -----------------
	160
	161	INTEGER i,j,l,ig0,ig,iq,iiq
	162	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
	163	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
	164	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
[1056]	165
[1549]	166	! REAL zrot(iip1,jjb_v:jje_v,llm) ! AdlC May 2014
	167	REAL :: zrot(iip1,jjm,llm)
[1]	168	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
	169	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
[8]	170	! ADAPTATION GCM POUR CP(T)
	171	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
	172	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
[1056]	173
[8]	174	! Ces calculs ne servent pas.
	175	! Si necessaire, decommenter ces variables et les calculs...
	176	! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
	177	! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
[1]	178	c
[1549]	179	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:), zrfi(:,:)
[1]	180	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
	181	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
	182	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
	183
	184	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
	185	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
[892]	186	REAL,ALLOCATABLE,SAVE :: zpk_omp(:,:)
[1]	187	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
	188	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
	189	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
	190	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
	191	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
[1549]	192	REAL,ALLOCATABLE,SAVE :: zrfi_omp(:,:)
[1]	193	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
	194	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
	195	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
	196	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
	197	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
	198	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
	199	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
	200	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
	201
	202	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	203	! Introduction du splitting (FH)
	204	! Question pour Yann :
	205	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
	206	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
	207	! soit allocatable (plutot par exemple que de passer une dimension
	208	! dépendant du process en argument des routines) et que, du coup,
	209	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
	210	! Tu confirmes ?
	211	! J'ai suivi le même principe pour les zdufic_omp
	212	! Mais c'est surement bien que tu controles.
	213	!
	214
	215	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
	216	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
	217	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
	218	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
	219	REAL jH_cur_split,zdt_split
	220	LOGICAL debut_split,lafin_split
	221	INTEGER isplit
	222	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	223
[892]	224	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp,
[1]	225	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
[1549]	226	c$OMP+ zrfi_omp,zqfi_omp,zdufi_omp,zdvfi_omp,
[1]	227	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
	228	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
	229
	230	LOGICAL,SAVE :: first_omp=.true.
	231	c$OMP THREADPRIVATE(first_omp)
	232
	233	REAL zsin(iim),zcos(iim),z1(iim)
	234	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	235	REAL unskap, pksurcp
[1056]	236	save unskap
	237
[1]	238	REAL SSUM
	239
[1189]	240	LOGICAL,SAVE :: firstcal=.true., debut=.true.
[1]	241	c$OMP THREADPRIVATE(firstcal,debut)
	242
	243	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
	244	INTEGER :: ierr
	245	#ifdef CPP_MPI
	246	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
	247	#else
	248	INTEGER,dimension(1,4) :: Status
	249	#endif
	250	INTEGER, dimension(4) :: Req
	251	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
	252	integer :: k,kstart,kend
	253	INTEGER :: offset
[1549]	254	INTEGER :: jjb,jje
[841]	255
[1056]	256	! For Titan only right now:
	257	! to allow for 2D computation of microphys and chemistry
	258	LOGICAL,save :: flag_moyzon
[1076]	259	REAL,allocatable,save :: tmpvar(:,:)
	260	REAL,allocatable,save :: tmpvarp1(:,:)
	261	REAL,allocatable,save :: tmpvarbar(:)
	262	REAL,allocatable,save :: tmpvarbarp1(:)
[1126]	263	real :: zz1,zz2
[1056]	264
[1]	265	c-----------------------------------------------------------------------
[1056]	266
[1]	267	c 1. Initialisations :
	268	c --------------------
	269
[1056]	270
[1]	271	klon=klon_mpi
	272
	273	IF ( firstcal ) THEN
	274	debut = .TRUE.
	275	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
	276	write(lunout,*) 'STOP dans calfis'
	277	write(lunout,*)
	278	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
	279	write(lunout,*) ' ngridmx jjm iim '
	280	write(lunout,*) ngridmx,jjm,iim
	281	STOP
	282	ENDIF
[1056]	283
	284	unskap = 1./ kappa
	285
[1076]	286	c$OMP MASTER
[1056]	287	flag_moyzon = .false.
	288	if(moyzon_ch.or.moyzon_mu) then
	289	flag_moyzon = .true.
[1071]	290	allocate(tmpvar(iip1,llm))
	291	allocate(tmpvarp1(iip1,llmp1))
	292	allocate(tmpvarbar(llm))
	293	allocate(tmpvarbarp1(llmp1))
[1056]	294	endif
	295
[1]	296	ALLOCATE(zpsrf(klon))
	297	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
	298	ALLOCATE(zphi(klon,llm),zphis(klon))
	299	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
	300	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
[8]	301	! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
	302	! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
[1549]	303	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm),zrfi(klon,llm))
[1]	304	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
	305	ALLOCATE(zdpsrf(klon))
	306	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
	307	ALLOCATE(flxwfi(klon,llm))
[8]	308	! ADAPTATION GCM POUR CP(T)
[52]	309	ALLOCATE(zteta(klon,llm))
[8]	310	ALLOCATE(zpk(klon,llm))
[1056]	311
[1251]	312	! zonal means. horizontal dimension should be iip1
[1256]	313	if (flag_moyzon) call moyzon_init(klon,llm,nqtot)
[1056]	314
	315	c------------------------------------------------------------------
	316	c moyennes globales pour les profils de pression et de temperature
[1126]	317	if(planet_type.eq."titan".or.planet_type.eq."venus") then
[1056]	318	call AllGather_Field(pp,iip1*jjp1,llmp1)
	319	call AllGather_Field(pteta,iip1*jjp1,llm)
	320	call AllGather_Field(ppk,iip1*jjp1,llm)
[1126]	321	call AllGather_Field(pphi,iip1*jjp1,llm)
	322	call AllGather_Field(pphis,iip1*jjp1,1)
[1056]	323	ALLOCATE(plevmoy(llm+1))
	324	ALLOCATE(playmoy(llm))
	325	ALLOCATE(tmoy(llm))
	326	ALLOCATE(tetamoy(llm))
	327	ALLOCATE(pkmoy(llm))
[1126]	328	ALLOCATE(phimoy(0:llm))
	329	ALLOCATE(zlevmoy(llm+1))
	330	ALLOCATE(zlaymoy(llm))
[1056]	331	plevmoy=0.
	332	do l=1,llmp1
	333	do i=1,iip1
	334	do j=1,jjp1
	335	plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
	336	enddo
	337	enddo
	338	enddo
	339	tetamoy=0.
	340	pkmoy=0.
[1126]	341	phimoy=0.
	342	do i=1,iip1
	343	do j=1,jjp1
	344	phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
	345	enddo
	346	enddo
[1056]	347	do l=1,llm
	348	do i=1,iip1
	349	do j=1,jjp1
	350	tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
	351	pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
[1126]	352	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
[1056]	353	enddo
	354	enddo
	355	enddo
[1126]	356	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
[1056]	357	call tpot2t_p(1,llm,tetamoy,tmoy,pkmoy)
[1126]	358	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
[1459]	359	zlaymoy(1:llm) = gradrad/(g*rad-phimoy(1:llm))-rad
[1126]	360	zlevmoy(1) = phimoy(0)/g
	361	DO l=2,llm
	362	zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
	363	zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l))
	364	zlevmoy(l)=(zz1zlaymoy(l-1)+zz2zlaymoy(l))/(zz1+zz2)
	365	ENDDO
	366	zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
[1071]	367	c-------------------
[1056]	368	c + lat index
[1071]	369	allocate(klat(klon))
	370	do ig0=1,klon
[1126]	371	j=index_j(ig0)
	372	klat(ig0)=j
[1071]	373	enddo
	374	endif ! planet_type=titan
[1056]	375	c------------------------------------------------------------------
	376
[1]	377	c$OMP END MASTER
	378	c$OMP BARRIER
	379	ELSE
	380	debut = .FALSE.
	381	ENDIF
	382
[1056]	383
[1]	384	c-----------------------------------------------------------------------
	385	c 40. transformation des variables dynamiques en variables physiques:
	386	c ---------------------------------------------------------------
	387
	388	c 41. pressions au sol (en Pascals)
	389	c ----------------------------------
	390
	391	c$OMP MASTER
	392	call start_timer(timer_physic)
	393	c$OMP END MASTER
	394
	395	c$OMP MASTER
	396	!CDIR ON_ADB(index_i)
	397	!CDIR ON_ADB(index_j)
	398	do ig0=1,klon
	399	i=index_i(ig0)
	400	j=index_j(ig0)
	401	zpsrf(ig0)=pps(i,j)
	402	enddo
	403	c$OMP END MASTER
	404
	405
[8]	406	c 42. pression intercouches et fonction d'Exner:
[1056]	407
[1]	408	c -----------------------------------------------------------------
	409	c .... zplev definis aux (llm +1) interfaces des couches ....
	410	c .... zplay definis aux ( llm ) milieux des couches ....
	411	c -----------------------------------------------------------------
	412
	413	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
[1056]	414
[1]	415	c print *,omp_rank,'klon--->',klon
	416	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	417	DO l = 1, llmp1
	418	!CDIR ON_ADB(index_i)
	419	!CDIR ON_ADB(index_j)
	420	do ig0=1,klon
	421	i=index_i(ig0)
	422	j=index_j(ig0)
	423	zplev( ig0,l ) = pp(i,j,l)
	424	enddo
	425	ENDDO
	426	c$OMP END DO NOWAIT
[1056]	427	if (flag_moyzon) then
	428	call AllGather_Field(pp,iip1*jjp1,llmp1)
	429	j=index_j(1)
	430	tmpvarp1(:,:) = pp(:,j,:)
	431	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
	432	zplevbar_mpi(1,:) = tmpvarbarp1
	433	do ig0=2,klon
	434	j=index_j(ig0)
	435	if (j.ne.index_j(ig0-1)) then
	436	tmpvarp1(:,:) = pp(:,j,:)
	437	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
	438	zplevbar_mpi(ig0,:) = tmpvarbarp1
	439	else
	440	zplevbar_mpi(ig0,:) = zplevbar_mpi(ig0-1,:)
	441	endif
	442	enddo
	443	endif
	444
[8]	445	! ADAPTATION GCM POUR CP(T)
[1056]	446	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[8]	447	DO l=1,llm
	448	!CDIR ON_ADB(index_i)
	449	!CDIR ON_ADB(index_j)
	450	do ig0=1,klon
	451	i=index_i(ig0)
	452	j=index_j(ig0)
	453	zpk(ig0,l)=ppk(i,j,l)
	454	zteta(ig0,l)=pteta(i,j,l)
	455	enddo
	456	ENDDO
	457	c$OMP END DO NOWAIT
[1056]	458	if (flag_moyzon) then
	459	call AllGather_Field(pteta,iip1*jjp1,llm)
	460	call AllGather_Field(ppk,iip1*jjp1,llm)
	461	j=index_j(1)
	462	tmpvar(:,:) = pteta(:,j,:)
	463	call moyzon(llm,tmpvar,tmpvarbar)
	464	ztetabar_mpi(1,:) = tmpvarbar
	465	tmpvar(:,:) = ppk(:,j,:)
	466	call moyzon(llm,tmpvar,tmpvarbar)
	467	zpkbar_mpi(1,:) = tmpvarbar
	468	call tpot2t_p(1,llm,ztetabar_mpi(1,:),ztfibar_mpi(1,:),
	469	& zpkbar_mpi(1,:))
	470	do ig0=2,klon
	471	j=index_j(ig0)
	472	if (j.ne.index_j(ig0-1)) then
	473	tmpvar(:,:) = pteta(:,j,:)
	474	call moyzon(llm,tmpvar,tmpvarbar)
	475	ztetabar_mpi(ig0,:) = tmpvarbar
	476	tmpvar(:,:) = ppk(:,j,:)
	477	call moyzon(llm,tmpvar,tmpvarbar)
	478	zpkbar_mpi(ig0,:) = tmpvarbar
	479	call tpot2t_p(1,llm,ztetabar_mpi(ig0,:),ztfibar_mpi(ig0,:),
	480	& zpkbar_mpi(ig0,:))
	481	else
	482	zpkbar_mpi(ig0,:) = zpkbar_mpi(ig0-1,:)
	483	ztetabar_mpi(ig0,:) = ztetabar_mpi(ig0-1,:)
	484	ztfibar_mpi(ig0,:) = ztfibar_mpi(ig0-1,:)
	485	endif
	486	enddo
	487	endif
[8]	488
[1]	489	c 43. temperature naturelle (en K) et pressions milieux couches .
	490	c ---------------------------------------------------------------
[8]	491
	492	! ADAPTATION GCM POUR CP(T)
[847]	493	call tpot2t_p(klon,llm,zteta,ztfi,zpk)
[8]	494
[1]	495	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	496	DO l=1,llm
	497	!CDIR ON_ADB(index_i)
	498	!CDIR ON_ADB(index_j)
	499	do ig0=1,klon
	500	i=index_i(ig0)
	501	j=index_j(ig0)
	502	pksurcp = ppk(i,j,l) / cpp
	503	zplay(ig0,l) = preff * pksurcp ** unskap
	504	enddo
	505	ENDDO
	506	c$OMP END DO NOWAIT
[1056]	507	if (flag_moyzon) then
	508	zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap
	509	endif
[1]	510
[8]	511	c 43.bis traceurs (tous intensifs)
[1]	512	c ---------------
	513
	514	DO iq=1,nqtot
	515	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	516	DO l=1,llm
	517	!CDIR ON_ADB(index_i)
	518	!CDIR ON_ADB(index_j)
	519	do ig0=1,klon
	520	i=index_i(ig0)
	521	j=index_j(ig0)
[1056]	522	zqfi(ig0,l,iq) = pq(i,j,l,iq)
[1]	523	enddo
	524	ENDDO
	525	c$OMP END DO NOWAIT
[8]	526	ENDDO ! of DO iq=1,nqtot
[1056]	527	if (flag_moyzon) then
	528	DO iq=1,nqtot
	529	call AllGather_Field(pq(:,:,:,iq),iip1*jjp1,llm)
	530	j=index_j(1)
	531	tmpvar(:,:) = pq(:,j,:,iq)
	532	call moyzon(llm,tmpvar,tmpvarbar)
	533	zqfibar_mpi(1,:,iq) = tmpvarbar
	534	do ig0=2,klon
	535	j=index_j(ig0)
	536	if (j.ne.index_j(ig0-1)) then
	537	tmpvar(:,:) = pq(:,j,:,iq)
	538	call moyzon(llm,tmpvar,tmpvarbar)
	539	zqfibar_mpi(ig0,:,iq) = tmpvarbar
	540	else
	541	zqfibar_mpi(ig0,:,iq) = zqfibar_mpi(ig0-1,:,iq)
	542	endif
	543	enddo
	544	ENDDO ! of DO iq=1,nqtot
	545	endif
[1]	546
	547
	548	c Geopotentiel calcule par rapport a la surface locale:
	549	c -----------------------------------------------------
	550
	551	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
	552
	553	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
	554
	555	c$OMP BARRIER
	556	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	557	DO l=1,llm
	558	DO ig=1,klon
	559	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	560	ENDDO
	561	ENDDO
	562	c$OMP END DO NOWAIT
	563
[1056]	564	if (flag_moyzon) then
	565	call AllGather_Field(pphis,iip1*jjp1,1)
	566	call AllGather_Field(pphi,iip1*jjp1,llm)
	567	j=index_j(1)
	568	tmpvar(:,1) = pphis(:,j)
	569	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
	570	zphisbar_mpi(1) = tmpvarbar(1)
	571	tmpvar(:,:) = pphi(:,j,:)
	572	call moyzon(llm,tmpvar,tmpvarbar)
	573	zphibar_mpi(1,:) = tmpvarbar-zphisbar_mpi(1)
	574	do ig0=2,klon
	575	j=index_j(ig0)
	576	if (j.ne.index_j(ig0-1)) then
	577	tmpvar(:,1) = pphis(:,j)
	578	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
	579	zphisbar_mpi(ig0) = tmpvarbar(1)
	580	tmpvar(:,:) = pphi(:,j,:)
	581	call moyzon(llm,tmpvar,tmpvarbar)
	582	zphibar_mpi(ig0,:) = tmpvarbar-zphisbar_mpi(ig0)
	583	else
	584	zphisbar_mpi(ig0) = zphisbar_mpi(ig0-1)
	585	zphibar_mpi(ig0,:) = zphibar_mpi(ig0-1,:)
	586	endif
	587	enddo
	588	endif
[1]	589
	590	c
	591	c 45. champ u:
	592	c ------------
	593
	594	kstart=1
	595	kend=klon
	596
[1572]	597	if (is_north_pole_dyn) kstart=2
	598	if (is_south_pole_dyn) kend=klon-1
[1]	599
	600	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	601	DO l=1,llm
	602	!CDIR ON_ADB(index_i)
	603	!CDIR ON_ADB(index_j)
	604	!CDIR SPARSE
	605	do ig0=kstart,kend
	606	i=index_i(ig0)
	607	j=index_j(ig0)
	608	if (i==1) then
	609	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
	610	$ + pucov(1,j,l)/cu(1,j) )
	611	else
	612	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
	613	$ + pucov(i,j,l)/cu(i,j) )
	614	endif
	615	enddo
	616	ENDDO
	617	c$OMP END DO NOWAIT
	618
[1549]	619	c
	620	C Alvaro de la Camara (May 2014)
	621	C 46.1 Calcul de la vorticite et passage sur la grille physique
	622	C --------------------------------------------------------------
	623
	624	jjb=jj_begin_dyn-1
	625	jje=jj_end_dyn+1
[1572]	626	if (is_north_pole_dyn) jjb=1
	627	if (is_south_pole_dyn) jje=jjm
[1549]	628
	629	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	630
	631	DO l=1,llm
	632	do i=1,iim
	633	do j=jjb,jje
	634	zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
	635	$ + pucov(i,j+1,l) - pucov(i,j,l))
	636	$ / (cu(i,j)+cu(i,j+1))
	637	$ / (cv(i+1,j)+cv(i,j)) *4
	638	enddo
	639	enddo
	640	ENDDO
	641
	642
	643	c 46.2champ v:
[1]	644	c -----------
	645
	646	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	647	DO l=1,llm
	648	!CDIR ON_ADB(index_i)
	649	!CDIR ON_ADB(index_j)
	650	DO ig0=kstart,kend
	651	i=index_i(ig0)
	652	j=index_j(ig0)
	653	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
	654	$ + pvcov(i,j,l)/cv(i,j) )
	655
[1549]	656	if (j==1 .OR. j==jjp1) then ! AdlC MAY 2014
	657	zrfi(ig0,l) = 0 ! AdlC MAY 2014
	658	else
	659	if(i==1)then
	660	zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
	661	$ +zrot(1,j-1,l)+zrot(1,j,l)) ! AdlC MAY 2014
	662	else
	663	zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
	664	$ +zrot(i,j-1,l)+zrot(i,j,l)) ! AdlC MAY 2014
	665	endif
	666	endif
	667
	668
[1]	669	ENDDO
	670	ENDDO
	671	c$OMP END DO NOWAIT
	672
	673	c 47. champs de vents aux pole nord
	674	c ------------------------------
	675	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	676	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	677
[1572]	678	if (is_north_pole_dyn) then
[1]	679	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	680	DO l=1,llm
	681
	682	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	683	DO i=2,iim
	684	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	685	ENDDO
	686
	687	DO i=1,iim
	688	zcos(i) = COS(rlonv(i))*z1(i)
	689	zsin(i) = SIN(rlonv(i))*z1(i)
	690	ENDDO
	691
	692	zufi(1,l) = SSUM(iim,zcos,1)/pi
	693	zvfi(1,l) = SSUM(iim,zsin,1)/pi
[1549]	694	zrfi(1,l) = 0.
[1]	695
	696	ENDDO
	697	c$OMP END DO NOWAIT
	698	endif
	699
	700
	701	c 48. champs de vents aux pole sud:
	702	c ---------------------------------
	703	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	704	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	705
[1572]	706	if (is_south_pole_dyn) then
[1]	707	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	708	DO l=1,llm
	709
	710	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	711	DO i=2,iim
	712	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	713	ENDDO
	714
	715	DO i=1,iim
	716	zcos(i) = COS(rlonv(i))*z1(i)
	717	zsin(i) = SIN(rlonv(i))*z1(i)
	718	ENDDO
	719
	720	zufi(klon,l) = SSUM(iim,zcos,1)/pi
	721	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
[1549]	722	zrfi(klon,l) = 0.
[1]	723	ENDDO
	724	c$OMP END DO NOWAIT
	725	endif
	726
	727	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
	728	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
	729
	730	c-----------------------------------------------------------------------
	731	c Appel de la physique:
	732	c ---------------------
	733
[8]	734	! Appel de la physique: pose probleme quand on tourne
	735	! SANS physique, car physiq.F est dans le repertoire phy[]...
	736	! Il faut une cle CPP_PHYS
[1]	737
[8]	738	! Le fait que les arguments de physiq soient differents selon les planetes
	739	! ne pose pas de probleme a priori.
	740
	741
[1]	742	c$OMP BARRIER
	743	if (first_omp) then
	744	klon=klon_omp
	745
	746	allocate(zplev_omp(klon,llm+1))
	747	allocate(zplay_omp(klon,llm))
[892]	748	allocate(zpk_omp(klon,llm))
[1]	749	allocate(zphi_omp(klon,llm))
	750	allocate(zphis_omp(klon))
	751	allocate(presnivs_omp(llm))
	752	allocate(zufi_omp(klon,llm))
	753	allocate(zvfi_omp(klon,llm))
[1549]	754	allocate(zrfi_omp(klon,llm)) ! LG Ari 2014
[1]	755	allocate(ztfi_omp(klon,llm))
	756	allocate(zqfi_omp(klon,llm,nqtot))
	757	allocate(zdufi_omp(klon,llm))
	758	allocate(zdvfi_omp(klon,llm))
	759	allocate(zdtfi_omp(klon,llm))
	760	allocate(zdqfi_omp(klon,llm,nqtot))
	761	allocate(zdufic_omp(klon,llm))
	762	allocate(zdvfic_omp(klon,llm))
	763	allocate(zdtfic_omp(klon,llm))
	764	allocate(zdqfic_omp(klon,llm,nqtot))
	765	allocate(zdpsrf_omp(klon))
	766	allocate(flxwfi_omp(klon,llm))
[1056]	767
[1251]	768	if (flag_moyzon) call moyzon_init_omp(klon,llm,nqtot)
[1056]	769
[1]	770	first_omp=.false.
	771	endif
	772
	773
	774	klon=klon_omp
	775	offset=klon_omp_begin-1
	776
	777	do l=1,llm+1
	778	do i=1,klon
	779	zplev_omp(i,l)=zplev(offset+i,l)
	780	enddo
	781	enddo
	782
	783	do l=1,llm
	784	do i=1,klon
	785	zplay_omp(i,l)=zplay(offset+i,l)
	786	enddo
	787	enddo
	788
[892]	789	do l=1,llm
	790	do i=1,klon
	791	zpk_omp(i,l)=zpk(offset+i,l)
	792	enddo
	793	enddo
	794
[1]	795	do l=1,llm
	796	do i=1,klon
	797	zphi_omp(i,l)=zphi(offset+i,l)
	798	enddo
	799	enddo
	800
	801	do i=1,klon
	802	zphis_omp(i)=zphis(offset+i)
	803	enddo
	804
	805
	806	do l=1,llm
	807	presnivs_omp(l)=presnivs(l)
	808	enddo
	809
	810	do l=1,llm
	811	do i=1,klon
	812	zufi_omp(i,l)=zufi(offset+i,l)
	813	enddo
	814	enddo
	815
	816	do l=1,llm
	817	do i=1,klon
	818	zvfi_omp(i,l)=zvfi(offset+i,l)
	819	enddo
	820	enddo
	821
	822	do l=1,llm
	823	do i=1,klon
[1549]	824	zrfi_omp(i,l)=zrfi(offset+i,l)
	825	enddo
	826	enddo
	827
	828
	829	do l=1,llm
	830	do i=1,klon
[1]	831	ztfi_omp(i,l)=ztfi(offset+i,l)
	832	enddo
	833	enddo
	834
	835	do iq=1,nqtot
	836	do l=1,llm
	837	do i=1,klon
	838	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
	839	enddo
	840	enddo
	841	enddo
	842
	843	do l=1,llm
	844	do i=1,klon
	845	zdufi_omp(i,l)=zdufi(offset+i,l)
	846	enddo
	847	enddo
	848
	849	do l=1,llm
	850	do i=1,klon
	851	zdvfi_omp(i,l)=zdvfi(offset+i,l)
	852	enddo
	853	enddo
	854
	855	do l=1,llm
	856	do i=1,klon
	857	zdtfi_omp(i,l)=zdtfi(offset+i,l)
	858	enddo
	859	enddo
	860
	861	do iq=1,nqtot
	862	do l=1,llm
	863	do i=1,klon
	864	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
	865	enddo
	866	enddo
	867	enddo
	868
	869	do i=1,klon
	870	zdpsrf_omp(i)=zdpsrf(offset+i)
	871	enddo
	872
	873	do l=1,llm
	874	do i=1,klon
	875	flxwfi_omp(i,l)=flxwfi(offset+i,l)
	876	enddo
	877	enddo
[1056]	878
	879	if (flag_moyzon) then
	880	do l=1,llm+1
	881	do i=1,klon
	882	zplevbar(i,l)=zplevbar_mpi(offset+i,l)
	883	enddo
	884	enddo
	885
	886	do l=1,llm
	887	do i=1,klon
	888	zplaybar(i,l)=zplaybar_mpi(offset+i,l)
	889	enddo
	890	enddo
	891
	892	do l=1,llm
	893	do i=1,klon
	894	zphibar(i,l)=zphibar_mpi(offset+i,l)
	895	enddo
	896	enddo
	897
	898	do i=1,klon
	899	zphisbar(i)=zphisbar_mpi(offset+i)
	900	enddo
	901
	902	do l=1,llm
	903	do i=1,klon
	904	ztfibar(i,l)=ztfibar_mpi(offset+i,l)
	905	enddo
	906	enddo
	907
	908	do iq=1,nqtot
	909	do l=1,llm
	910	do i=1,klon
	911	zqfibar(i,l,iq)=zqfibar_mpi(offset+i,l,iq)
	912	enddo
	913	enddo
	914	enddo
	915	endif
	916
[1]	917
	918	c$OMP BARRIER
	919
	920	!$OMP MASTER
	921	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
	922	!$OMP END MASTER
	923	zdt_split=dtphys/nsplit_phys
	924	zdufic_omp(:,:)=0.
	925	zdvfic_omp(:,:)=0.
	926	zdtfic_omp(:,:)=0.
	927	zdqfic_omp(:,:,:)=0.
	928
[776]	929	#ifdef CPP_PHYS
[1]	930	do isplit=1,nsplit_phys
	931
	932	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
	933	debut_split=debut.and.isplit==1
	934	lafin_split=lafin.and.isplit==nsplit_phys
	935
[1549]	936	CALL call_physiq(klon,llm,nqtot,tname,
	937	& debut_split,lafin_split,
	938	& jD_cur,jH_cur_split,zdt_split,
	939	& zplev_omp,zplay_omp,
	940	& zpk_omp,zphi_omp,zphis_omp,
	941	& presnivs_omp,
	942	& zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
	943	& flxwfi_omp,pducov,
	944	& zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
[1576]	945	& zdpsrf_omp)
[1]	946
[1549]	947	! if (planet_type=="earth") then
	948	! CALL physiq (klon,
	949	! . llm,
	950	! . debut_split,
	951	! . lafin_split,
	952	! . jD_cur,
	953	! . jH_cur_split,
	954	! . zdt_split,
	955	! . zplev_omp,
	956	! . zplay_omp,
	957	! . zphi_omp,
	958	! . zphis_omp,
	959	! . presnivs_omp,
	960	! . zufi_omp,
	961	! . zvfi_omp,
	962	! . ztfi_omp,
	963	! . zqfi_omp,
	964	! . flxwfi_omp,
	965	! . zdufi_omp,
	966	! . zdvfi_omp,
	967	! . zdtfi_omp,
	968	! . zdqfi_omp,
	969	! . zdpsrf_omp,
	970	! . pducov)
	971	!
	972	! else if ( planet_type=="generic" ) then
	973	!
	974	! CALL physiq (klon, !! ngrid
	975	! . llm, !! nlayer
	976	! . nqtot, !! nq
	977	! . tname, !! tracer names from dynamical core (given in infotrac)
	978	! . debut_split, !! firstcall
	979	! . lafin_split, !! lastcall
	980	! . jD_cur, !! pday. see leapfrog_p
	981	! . jH_cur_split, !! ptime "fraction of day"
	982	! . zdt_split, !! ptimestep
	983	! . zplev_omp, !! pplev
	984	! . zplay_omp, !! pplay
	985	! . zphi_omp, !! pphi
	986	! . zufi_omp, !! pu
	987	! . zvfi_omp, !! pv
	988	! . ztfi_omp, !! pt
	989	! . zqfi_omp, !! pq
	990	! . flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
	991	! . zdufi_omp, !! pdu
	992	! . zdvfi_omp, !! pdv
	993	! . zdtfi_omp, !! pdt
	994	! . zdqfi_omp, !! pdq
	995	! . zdpsrf_omp, !! pdpsrf
	996	! . tracerdyn) !! tracerdyn <-- utilite ???
	997	!
	998	! else if ( planet_type=="mars" ) then
	999	!
	1000	! CALL physiq (klon, ! ngrid
	1001	! . llm, ! nlayer
	1002	! . nqtot, ! nq
	1003	! . debut_split, ! firstcall
	1004	! . lafin_split, ! lastcall
	1005	! . jD_cur, ! pday
	1006	! . jH_cur_split, ! ptime
	1007	! . zdt_split, ! ptimestep
	1008	! . zplev_omp, ! pplev
	1009	! . zplay_omp, ! pplay
	1010	! . zphi_omp, ! pphi
	1011	! . zufi_omp, ! pu
	1012	! . zvfi_omp, ! pv
	1013	! . ztfi_omp, ! pt
	1014	! . zqfi_omp, ! pq
	1015	! . flxwfi_omp, ! pw
	1016	! . zdufi_omp, ! pdu
	1017	! . zdvfi_omp, ! pdv
	1018	! . zdtfi_omp, ! pdt
	1019	! . zdqfi_omp, ! pdq
	1020	! . zdpsrf_omp, ! pdpsrf
	1021	! . tracerdyn) ! tracerdyn (somewhat obsolete)
	1022	!
	1023	! else if ((planet_type=="titan").or.(planet_type=="venus")) then
	1024	!
	1025	! CALL physiq (klon,
	1026	! . llm,
	1027	! . nqtot,
	1028	! . debut_split,
	1029	! . lafin_split,
	1030	! . jD_cur,
	1031	! . jH_cur_split,
	1032	! . zdt_split,
	1033	! . zplev_omp,
	1034	! . zplay_omp,
	1035	! . zpk_omp,
	1036	! . zphi_omp,
	1037	! . zphis_omp,
	1038	! . presnivs_omp,
	1039	! . zufi_omp,
	1040	! . zvfi_omp,
	1041	! . ztfi_omp,
	1042	! . zqfi_omp,
	1043	! . flxwfi_omp,
	1044	! . zdufi_omp,
	1045	! . zdvfi_omp,
	1046	! . zdtfi_omp,
	1047	! . zdqfi_omp,
	1048	! . zdpsrf_omp)
	1049	!
	1050	! else ! unknown "planet_type"
	1051	!
	1052	! write(lunout,*) "calfis_p: error, unknown planet_type: ",
	1053	! & trim(planet_type)
	1054	! stop
	1055	!
	1056	! endif ! planet_type
[1]	1057	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
	1058	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
	1059	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
	1060	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
	1061
	1062	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
	1063	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
	1064	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
	1065	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
	1066
	1067	enddo
	1068
[776]	1069	#endif
	1070	! of #ifdef CPP_PHYS
	1071
[1]	1072	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
	1073	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
	1074	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
	1075	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
	1076
	1077	c$OMP BARRIER
	1078
	1079	do l=1,llm+1
	1080	do i=1,klon
	1081	zplev(offset+i,l)=zplev_omp(i,l)
	1082	enddo
	1083	enddo
	1084
	1085	do l=1,llm
	1086	do i=1,klon
	1087	zplay(offset+i,l)=zplay_omp(i,l)
	1088	enddo
	1089	enddo
	1090
	1091	do l=1,llm
	1092	do i=1,klon
	1093	zphi(offset+i,l)=zphi_omp(i,l)
	1094	enddo
	1095	enddo
	1096
	1097
	1098	do i=1,klon
	1099	zphis(offset+i)=zphis_omp(i)
	1100	enddo
	1101
	1102
	1103	do l=1,llm
	1104	presnivs(l)=presnivs_omp(l)
	1105	enddo
	1106
	1107	do l=1,llm
	1108	do i=1,klon
	1109	zufi(offset+i,l)=zufi_omp(i,l)
	1110	enddo
	1111	enddo
	1112
	1113	do l=1,llm
	1114	do i=1,klon
	1115	zvfi(offset+i,l)=zvfi_omp(i,l)
	1116	enddo
	1117	enddo
	1118
	1119	do l=1,llm
	1120	do i=1,klon
	1121	ztfi(offset+i,l)=ztfi_omp(i,l)
	1122	enddo
	1123	enddo
	1124
	1125	do iq=1,nqtot
	1126	do l=1,llm
	1127	do i=1,klon
	1128	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
	1129	enddo
	1130	enddo
	1131	enddo
	1132
	1133	do l=1,llm
	1134	do i=1,klon
	1135	zdufi(offset+i,l)=zdufi_omp(i,l)
	1136	enddo
	1137	enddo
	1138
	1139	do l=1,llm
	1140	do i=1,klon
	1141	zdvfi(offset+i,l)=zdvfi_omp(i,l)
	1142	enddo
	1143	enddo
	1144
	1145	do l=1,llm
	1146	do i=1,klon
	1147	zdtfi(offset+i,l)=zdtfi_omp(i,l)
	1148	enddo
	1149	enddo
	1150
	1151	do iq=1,nqtot
	1152	do l=1,llm
	1153	do i=1,klon
	1154	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
	1155	enddo
	1156	enddo
	1157	enddo
	1158
	1159	do i=1,klon
	1160	zdpsrf(offset+i)=zdpsrf_omp(i)
	1161	enddo
	1162
	1163
	1164	klon=klon_mpi
	1165	500 CONTINUE
	1166	c$OMP BARRIER
	1167
	1168	c$OMP MASTER
	1169	call stop_timer(timer_physic)
	1170	c$OMP END MASTER
	1171
	1172	IF (using_mpi) THEN
	1173
	1174	if (MPI_rank>0) then
	1175
	1176	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1177	DO l=1,llm
	1178	du_send(1:iim,l)=zdufi(1:iim,l)
	1179	dv_send(1:iim,l)=zdvfi(1:iim,l)
	1180	ENDDO
	1181	c$OMP END DO NOWAIT
	1182
	1183	c$OMP BARRIER
	1184	#ifdef CPP_MPI
	1185	c$OMP MASTER
	1186	!$OMP CRITICAL (MPI)
	1187	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
	1188	& COMM_LMDZ,Req(1),ierr)
	1189	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
	1190	& COMM_LMDZ,Req(2),ierr)
	1191	!$OMP END CRITICAL (MPI)
	1192	c$OMP END MASTER
	1193	#endif
	1194	c$OMP BARRIER
	1195
	1196	endif
	1197
	1198	if (MPI_rank<MPI_Size-1) then
	1199	c$OMP BARRIER
	1200	#ifdef CPP_MPI
	1201	c$OMP MASTER
	1202	!$OMP CRITICAL (MPI)
	1203	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
	1204	& COMM_LMDZ,Req(3),ierr)
	1205	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
	1206	& COMM_LMDZ,Req(4),ierr)
	1207	!$OMP END CRITICAL (MPI)
	1208	c$OMP END MASTER
	1209	#endif
	1210	endif
	1211
	1212	c$OMP BARRIER
	1213
	1214
	1215	#ifdef CPP_MPI
	1216	c$OMP MASTER
	1217	!$OMP CRITICAL (MPI)
	1218	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
	1219	call MPI_WAITALL(4,Req(1),Status,ierr)
	1220	else if (MPI_rank>0) then
	1221	call MPI_WAITALL(2,Req(1),Status,ierr)
	1222	else if (MPI_rank <MPI_Size-1) then
	1223	call MPI_WAITALL(2,Req(3),Status,ierr)
	1224	endif
	1225	!$OMP END CRITICAL (MPI)
	1226	c$OMP END MASTER
	1227	#endif
	1228
	1229	c$OMP BARRIER
	1230
	1231	ENDIF ! using_mpi
	1232
	1233
	1234	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1235	DO l=1,llm
	1236
	1237	zdufi2(1:klon,l)=zdufi(1:klon,l)
	1238	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
	1239
	1240	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
	1241	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
	1242
	1243	pdhfi(:,jj_begin,l)=0
	1244	pdqfi(:,jj_begin,l,:)=0
	1245	pdufi(:,jj_begin,l)=0
	1246	pdvfi(:,jj_begin,l)=0
	1247
[1572]	1248	if (.not. is_south_pole_dyn) then
[1]	1249	pdhfi(:,jj_end,l)=0
	1250	pdqfi(:,jj_end,l,:)=0
	1251	pdufi(:,jj_end,l)=0
	1252	pdvfi(:,jj_end,l)=0
	1253	endif
	1254
	1255	ENDDO
	1256	c$OMP END DO NOWAIT
	1257
	1258	c$OMP MASTER
	1259	pdpsfi(:,jj_begin)=0
[1572]	1260	if (.not. is_south_pole_dyn) then
[1]	1261	pdpsfi(:,jj_end)=0
	1262	endif
	1263	c$OMP END MASTER
	1264	c-----------------------------------------------------------------------
	1265	c transformation des tendances physiques en tendances dynamiques:
	1266	c ---------------------------------------------------------------
	1267
	1268	c tendance sur la pression :
	1269	c -----------------------------------
	1270	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
	1271	c
	1272	c 62. enthalpie potentielle
	1273	c ---------------------
[8]	1274
[1]	1275	kstart=1
	1276	kend=klon
	1277
[1572]	1278	if (is_north_pole_dyn) kstart=2
	1279	if (is_south_pole_dyn) kend=klon-1
[8]	1280
	1281	! ADAPTATION GCM POUR CP(T)
[847]	1282	call t2tpot_p(klon,llm,ztfi,zteta,zpk)
[1]	1283
[8]	1284
[1]	1285	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1286	DO l=1,llm
	1287	!CDIR ON_ADB(index_i)
	1288	!CDIR ON_ADB(index_j)
	1289	!cdir NODEP
	1290	do ig0=kstart,kend
	1291	i=index_i(ig0)
	1292	j=index_j(ig0)
[8]	1293	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	1294	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
	1295	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	1296	if (i==1) then
	1297	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
	1298	endif
[1]	1299	enddo
	1300
[1572]	1301	if (is_north_pole_dyn) then
[1]	1302	DO i=1,iip1
[8]	1303	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	1304	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
[1]	1305	enddo
	1306	endif
	1307
[1572]	1308	if (is_south_pole_dyn) then
[1]	1309	DO i=1,iip1
[8]	1310	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
	1311	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
[1]	1312	ENDDO
	1313	endif
	1314	ENDDO
	1315	c$OMP END DO NOWAIT
	1316
	1317	c 62. humidite specifique
	1318	c ---------------------
	1319	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	1320	! DO iq=1,nqtot
	1321	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1322	! DO l=1,llm
	1323	!!!cdir NODEP
	1324	! do ig0=kstart,kend
	1325	! i=index_i(ig0)
	1326	! j=index_j(ig0)
	1327	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
	1328	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
	1329	! enddo
	1330	!
[1572]	1331	! if (is_north_pole_dyn) then
[1]	1332	! do i=1,iip1
	1333	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	1334	! enddo
	1335	! endif
	1336	!
[1572]	1337	! if (is_south_pole_dyn) then
[1]	1338	! do i=1,iip1
	1339	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
	1340	! enddo
	1341	! endif
	1342	! ENDDO
	1343	!c$OMP END DO NOWAIT
	1344	! ENDDO
	1345
[8]	1346	c 63. traceurs (tous en intensifs)
[1]	1347	c ------------
	1348	C initialisation des tendances
	1349
	1350	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1351	DO l=1,llm
	1352	pdqfi(:,:,l,:)=0.
	1353	ENDDO
	1354	c$OMP END DO NOWAIT
	1355
	1356	C
	1357	!cdir NODEP
	1358	DO iq=1,nqtot
	1359	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1360	DO l=1,llm
	1361	!CDIR ON_ADB(index_i)
	1362	!CDIR ON_ADB(index_j)
	1363	!cdir NODEP
	1364	DO ig0=kstart,kend
	1365	i=index_i(ig0)
	1366	j=index_j(ig0)
[1056]	1367	pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
	1368	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
[1]	1369	ENDDO
	1370
[1572]	1371	IF (is_north_pole_dyn) then
[1]	1372	DO i=1,iip1
[1056]	1373	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
[1]	1374	ENDDO
	1375	ENDIF
	1376
[1572]	1377	IF (is_south_pole_dyn) then
[1]	1378	DO i=1,iip1
[1056]	1379	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
[1]	1380	ENDDO
	1381	ENDIF
	1382
	1383	ENDDO
	1384	c$OMP END DO NOWAIT
	1385	ENDDO
	1386
	1387	c 65. champ u:
	1388	c ------------
	1389	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1390	DO l=1,llm
	1391	!CDIR ON_ADB(index_i)
	1392	!CDIR ON_ADB(index_j)
	1393	!cdir NODEP
	1394	do ig0=kstart,kend
	1395	i=index_i(ig0)
	1396	j=index_j(ig0)
	1397
	1398	if (i/=iim) then
	1399	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
	1400	endif
	1401
	1402	if (i==1) then
	1403	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
	1404	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
	1405	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
	1406	endif
	1407
	1408	enddo
	1409
[1572]	1410	if (is_north_pole_dyn) then
[1]	1411	DO i=1,iip1
	1412	pdufi(i,1,l) = 0.
	1413	ENDDO
	1414	endif
	1415
[1572]	1416	if (is_south_pole_dyn) then
[1]	1417	DO i=1,iip1
	1418	pdufi(i,jjp1,l) = 0.
	1419	ENDDO
	1420	endif
	1421
	1422	ENDDO
	1423	c$OMP END DO NOWAIT
	1424
	1425	c 67. champ v:
	1426	c ------------
	1427
	1428	kstart=1
	1429	kend=klon
	1430
[1572]	1431	if (is_north_pole_dyn) kstart=2
	1432	if (is_south_pole_dyn) kend=klon-1-iim
[1]	1433
	1434	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1435	DO l=1,llm
	1436	!CDIR ON_ADB(index_i)
	1437	!CDIR ON_ADB(index_j)
	1438	!cdir NODEP
	1439	do ig0=kstart,kend
	1440	i=index_i(ig0)
	1441	j=index_j(ig0)
	1442	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
	1443	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
	1444	$ zdvfi2(ig0+iim,l))
	1445	$ *cv(i,j)
	1446	enddo
	1447
	1448	ENDDO
	1449	c$OMP END DO NOWAIT
	1450
	1451
	1452	c 68. champ v pres des poles:
	1453	c ---------------------------
	1454	c v = U * cos(long) + V * SIN(long)
	1455
[1572]	1456	if (is_north_pole_dyn) then
[1]	1457
	1458	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1459	DO l=1,llm
	1460
	1461	DO i=1,iim
	1462	pdvfi(i,1,l)=
	1463	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	1464
	1465	pdvfi(i,1,l)=
	1466	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	1467	ENDDO
	1468
	1469	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	1470
	1471	ENDDO
	1472	c$OMP END DO NOWAIT
	1473
	1474	endif
	1475
[1572]	1476	if (is_south_pole_dyn) then
[1]	1477
	1478	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1479	DO l=1,llm
	1480
	1481	DO i=1,iim
	1482	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
	1483	$ +zdvfi(klon,l)*SIN(rlonv(i))
	1484
	1485	pdvfi(i,jjm,l)=
	1486	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
	1487	ENDDO
	1488
	1489	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	1490
	1491	ENDDO
	1492	c$OMP END DO NOWAIT
	1493
	1494	endif
	1495	c-----------------------------------------------------------------------
	1496
	1497	700 CONTINUE
	1498
	1499	firstcal = .FALSE.
	1500
	1501	#else
	1502	write(lunout,*)
	1503	& "calfis_p: for now can only work with parallel physics"
	1504	stop
	1505	#endif
[8]	1506	! of #ifdef CPP_PHYS
[1403]	1507	#endif
	1508	! of #ifdef CPP_PARA
[1]	1509	END

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

Download in other formats:

Original Format