Context Navigation

calfis_p.F @ 1243

Last change on this file since 1243 was 1189, checked in by emillour, 11 years ago

Common dynamics: a couple of bug fixes:

Correctly account for the change in pressure, mass, etc. after modifying surface pressure following a call to the physics.
Corrected tracer advection, which is computed using values at the beginning of the time step, so it is done at Matsuno forward step.

File size: 40.4 KB

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

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

Download in other formats:

Original Format