Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

calfis_p.F

Last change on this file 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

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_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
34	USE mod_interface_dyn_phys
35	! USE IOPHY
36	#endif
37	#ifdef CPP_PARA
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
41	USE Write_Field
42	Use Write_field_p
43	USE Times
44	USE cpdet_mod, only: tpot2t_p, t2tpot_p
45	! used only for zonal averages
46	USE moyzon_mod
47	#endif
48	USE infotrac, ONLY: nqtot, niadv, tname
49	USE control_mod, ONLY: planet_type, nsplit_phys
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
53	#ifdef CPP_PHYS
54	USE callphysiq_mod, ONLY: call_physiq
55	#endif
56
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
99	c pdhfi tendency for the potential temperature (K/s)
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
112	include "dimensions.h"
113	include "paramet.h"
114
115	INTEGER ngridmx
116	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
117
118	include "comgeom2.h"
119	include "iniprint.h"
120	#ifdef CPP_MPI
121	include 'mpif.h'
122	#endif
123	c Arguments :
124	c -----------
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
134
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
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....
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...
143
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
147	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
148
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)
155
156	#ifdef CPP_PARA
157	#ifdef CPP_PHYS
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(:)
165
166	! REAL zrot(iip1,jjb_v:jje_v,llm) ! AdlC May 2014
167	REAL :: zrot(iip1,jjm,llm)
168	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
169	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
170	! ADAPTATION GCM POUR CP(T)
171	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
172	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
173
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(:,:,:)
178	c
179	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:), zrfi(:,:)
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(:,:)
186	REAL,ALLOCATABLE,SAVE :: zpk_omp(:,:)
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(:,:)
192	REAL,ALLOCATABLE,SAVE :: zrfi_omp(:,:)
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
224	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp,
225	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
226	c$OMP+ zrfi_omp,zqfi_omp,zdufi_omp,zdvfi_omp,
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
236	save unskap
237
238	REAL SSUM
239
240	LOGICAL,SAVE :: firstcal=.true., debut=.true.
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
254	INTEGER :: jjb,jje
255
256	! For Titan only right now:
257	! to allow for 2D computation of microphys and chemistry
258	LOGICAL,save :: flag_moyzon
259	REAL,allocatable,save :: tmpvar(:,:)
260	REAL,allocatable,save :: tmpvarp1(:,:)
261	REAL,allocatable,save :: tmpvarbar(:)
262	REAL,allocatable,save :: tmpvarbarp1(:)
263	real :: zz1,zz2
264
265	c-----------------------------------------------------------------------
266
267	c 1. Initialisations :
268	c --------------------
269
270
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
283
284	unskap = 1./ kappa
285
286	c$OMP MASTER
287	flag_moyzon = .false.
288	if(moyzon_ch.or.moyzon_mu) then
289	flag_moyzon = .true.
290	allocate(tmpvar(iip1,llm))
291	allocate(tmpvarp1(iip1,llmp1))
292	allocate(tmpvarbar(llm))
293	allocate(tmpvarbarp1(llmp1))
294	endif
295
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))
301	! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
302	! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
303	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm),zrfi(klon,llm))
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))
308	! ADAPTATION GCM POUR CP(T)
309	ALLOCATE(zteta(klon,llm))
310	ALLOCATE(zpk(klon,llm))
311
312	! zonal means. horizontal dimension should be iip1
313	if (flag_moyzon) call moyzon_init(klon,llm,nqtot)
314
315	c------------------------------------------------------------------
316	c moyennes globales pour les profils de pression et de temperature
317	if(planet_type.eq."titan".or.planet_type.eq."venus") then
318	call AllGather_Field(pp,iip1*jjp1,llmp1)
319	call AllGather_Field(pteta,iip1*jjp1,llm)
320	call AllGather_Field(ppk,iip1*jjp1,llm)
321	call AllGather_Field(pphi,iip1*jjp1,llm)
322	call AllGather_Field(pphis,iip1*jjp1,1)
323	ALLOCATE(plevmoy(llm+1))
324	ALLOCATE(playmoy(llm))
325	ALLOCATE(tmoy(llm))
326	ALLOCATE(tetamoy(llm))
327	ALLOCATE(pkmoy(llm))
328	ALLOCATE(phimoy(0:llm))
329	ALLOCATE(zlevmoy(llm+1))
330	ALLOCATE(zlaymoy(llm))
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.
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
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)
352	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
353	enddo
354	enddo
355	enddo
356	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
357	call tpot2t_p(1,llm,tetamoy,tmoy,pkmoy)
358	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
359	zlaymoy(1:llm) = gradrad/(g*rad-phimoy(1:llm))-rad
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))
367	c-------------------
368	c + lat index
369	allocate(klat(klon))
370	do ig0=1,klon
371	j=index_j(ig0)
372	klat(ig0)=j
373	enddo
374	endif ! planet_type=titan
375	c------------------------------------------------------------------
376
377	c$OMP END MASTER
378	c$OMP BARRIER
379	ELSE
380	debut = .FALSE.
381	ENDIF
382
383
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
406	c 42. pression intercouches et fonction d'Exner:
407
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 ....
414
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
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
445	! ADAPTATION GCM POUR CP(T)
446	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
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
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
488
489	c 43. temperature naturelle (en K) et pressions milieux couches .
490	c ---------------------------------------------------------------
491
492	! ADAPTATION GCM POUR CP(T)
493	call tpot2t_p(klon,llm,zteta,ztfi,zpk)
494
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
507	if (flag_moyzon) then
508	zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap
509	endif
510
511	c 43.bis traceurs (tous intensifs)
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)
522	zqfi(ig0,l,iq) = pq(i,j,l,iq)
523	enddo
524	ENDDO
525	c$OMP END DO NOWAIT
526	ENDDO ! of DO iq=1,nqtot
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
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
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
589
590	c
591	c 45. champ u:
592	c ------------
593
594	kstart=1
595	kend=klon
596
597	if (is_north_pole_dyn) kstart=2
598	if (is_south_pole_dyn) kend=klon-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
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
626	if (is_north_pole_dyn) jjb=1
627	if (is_south_pole_dyn) jje=jjm
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:
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
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
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
678	if (is_north_pole_dyn) then
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
694	zrfi(1,l) = 0.
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
706	if (is_south_pole_dyn) then
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
722	zrfi(klon,l) = 0.
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
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
737
738	! Le fait que les arguments de physiq soient differents selon les planetes
739	! ne pose pas de probleme a priori.
740
741
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))
748	allocate(zpk_omp(klon,llm))
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))
754	allocate(zrfi_omp(klon,llm)) ! LG Ari 2014
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))
767
768	if (flag_moyzon) call moyzon_init_omp(klon,llm,nqtot)
769
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
789	do l=1,llm
790	do i=1,klon
791	zpk_omp(i,l)=zpk(offset+i,l)
792	enddo
793	enddo
794
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
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
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
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
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
929	#ifdef CPP_PHYS
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
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,
945	& zdpsrf_omp)
946
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
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
1069	#endif
1070	! of #ifdef CPP_PHYS
1071
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
1248	if (.not. is_south_pole_dyn) then
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
1260	if (.not. is_south_pole_dyn) then
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 ---------------------
1274
1275	kstart=1
1276	kend=klon
1277
1278	if (is_north_pole_dyn) kstart=2
1279	if (is_south_pole_dyn) kend=klon-1
1280
1281	! ADAPTATION GCM POUR CP(T)
1282	call t2tpot_p(klon,llm,ztfi,zteta,zpk)
1283
1284
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)
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
1299	enddo
1300
1301	if (is_north_pole_dyn) then
1302	DO i=1,iip1
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
1305	enddo
1306	endif
1307
1308	if (is_south_pole_dyn) then
1309	DO i=1,iip1
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
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	!
1331	! if (is_north_pole_dyn) then
1332	! do i=1,iip1
1333	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1334	! enddo
1335	! endif
1336	!
1337	! if (is_south_pole_dyn) then
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
1346	c 63. traceurs (tous en intensifs)
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)
1367	pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1368	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1369	ENDDO
1370
1371	IF (is_north_pole_dyn) then
1372	DO i=1,iip1
1373	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1374	ENDDO
1375	ENDIF
1376
1377	IF (is_south_pole_dyn) then
1378	DO i=1,iip1
1379	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
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
1410	if (is_north_pole_dyn) then
1411	DO i=1,iip1
1412	pdufi(i,1,l) = 0.
1413	ENDDO
1414	endif
1415
1416	if (is_south_pole_dyn) then
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
1431	if (is_north_pole_dyn) kstart=2
1432	if (is_south_pole_dyn) kend=klon-1-iim
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
1456	if (is_north_pole_dyn) then
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
1476	if (is_south_pole_dyn) then
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
1506	! of #ifdef CPP_PHYS
1507	#endif
1508	! of #ifdef CPP_PARA
1509	END

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

Download in other formats:

Original Format