Context Navigation

calfis_p.F @ 1572

Last change on this file since 1572 was 1572, checked in by emillour, 8 years ago

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2500 of LMDZ5)

arch:
remove ifort debug option '-check all', replace it with '-check bounds,format,output_conversion,pointers,uninit' (i.e. get it to stop complaining about copying into temporary arrays)

dyn3d_common:
comconst_mod.F90 : add ngroup

dyn3d:
gcm.F90 : minor bug fix (arguments to a call_abort())
leapfrog.F90 : recompute geopotential for bilan_dyn outputs
conf_gcm.F90 : read "ngroup" from run.def
groupe.F , groupeun.F : ngroup no longer a local parameter

dyn3d_par:
conf_gcm.F90 : read "ngroup" from run.def
groupe_p.F , groupeun_p.F : ngroup no longer a local parameter

misc:
regr1_step_av_m.F90 : removed (not used)

phy_common:
mod_phys_lmdz_mpi_transfert.F90 , mod_phys_lmdz_mpi_data.F90 : change is_north_pole and is_south_pole to is_north_pole_dyn and is_south_pole_dyn
mod_phys_lmdz_omp_data.F90 : introduce is_nort_pole_phy and is_south_pole_phy

dynphy_lonlat:
mod_interface_dyn_phys.F90 : use is_north_pole_dyn and is_south_pole_dyn
calfis_p.F : use is_north_pole_dyn and is_south_pole_dyn

phyvenus:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

phytitan:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

File size: 42.5 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	LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
257
258	! For Titan only right now:
259	! to allow for 2D computation of microphys and chemistry
260	LOGICAL,save :: flag_moyzon
261	REAL,allocatable,save :: tmpvar(:,:)
262	REAL,allocatable,save :: tmpvarp1(:,:)
263	REAL,allocatable,save :: tmpvarbar(:)
264	REAL,allocatable,save :: tmpvarbarp1(:)
265	real :: zz1,zz2
266
267	c-----------------------------------------------------------------------
268
269	c 1. Initialisations :
270	c --------------------
271
272
273	klon=klon_mpi
274
275	IF ( firstcal ) THEN
276	debut = .TRUE.
277	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
278	write(lunout,*) 'STOP dans calfis'
279	write(lunout,*)
280	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
281	write(lunout,*) ' ngridmx jjm iim '
282	write(lunout,*) ngridmx,jjm,iim
283	STOP
284	ENDIF
285
286	unskap = 1./ kappa
287
288	c$OMP MASTER
289	flag_moyzon = .false.
290	if(moyzon_ch.or.moyzon_mu) then
291	flag_moyzon = .true.
292	allocate(tmpvar(iip1,llm))
293	allocate(tmpvarp1(iip1,llmp1))
294	allocate(tmpvarbar(llm))
295	allocate(tmpvarbarp1(llmp1))
296	endif
297
298	ALLOCATE(zpsrf(klon))
299	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
300	ALLOCATE(zphi(klon,llm),zphis(klon))
301	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
302	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
303	! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
304	! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
305	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm),zrfi(klon,llm))
306	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
307	ALLOCATE(zdpsrf(klon))
308	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
309	ALLOCATE(flxwfi(klon,llm))
310	! ADAPTATION GCM POUR CP(T)
311	ALLOCATE(zteta(klon,llm))
312	ALLOCATE(zpk(klon,llm))
313
314	! zonal means. horizontal dimension should be iip1
315	if (flag_moyzon) call moyzon_init(klon,llm,nqtot)
316
317	c------------------------------------------------------------------
318	c moyennes globales pour les profils de pression et de temperature
319	if(planet_type.eq."titan".or.planet_type.eq."venus") then
320	call AllGather_Field(pp,iip1*jjp1,llmp1)
321	call AllGather_Field(pteta,iip1*jjp1,llm)
322	call AllGather_Field(ppk,iip1*jjp1,llm)
323	call AllGather_Field(pphi,iip1*jjp1,llm)
324	call AllGather_Field(pphis,iip1*jjp1,1)
325	ALLOCATE(plevmoy(llm+1))
326	ALLOCATE(playmoy(llm))
327	ALLOCATE(tmoy(llm))
328	ALLOCATE(tetamoy(llm))
329	ALLOCATE(pkmoy(llm))
330	ALLOCATE(phimoy(0:llm))
331	ALLOCATE(zlevmoy(llm+1))
332	ALLOCATE(zlaymoy(llm))
333	plevmoy=0.
334	do l=1,llmp1
335	do i=1,iip1
336	do j=1,jjp1
337	plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
338	enddo
339	enddo
340	enddo
341	tetamoy=0.
342	pkmoy=0.
343	phimoy=0.
344	do i=1,iip1
345	do j=1,jjp1
346	phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
347	enddo
348	enddo
349	do l=1,llm
350	do i=1,iip1
351	do j=1,jjp1
352	tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
353	pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
354	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
355	enddo
356	enddo
357	enddo
358	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
359	call tpot2t_p(1,llm,tetamoy,tmoy,pkmoy)
360	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
361	zlaymoy(1:llm) = gradrad/(g*rad-phimoy(1:llm))-rad
362	zlevmoy(1) = phimoy(0)/g
363	DO l=2,llm
364	zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
365	zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l))
366	zlevmoy(l)=(zz1zlaymoy(l-1)+zz2zlaymoy(l))/(zz1+zz2)
367	ENDDO
368	zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
369	c-------------------
370	c + lat index
371	allocate(klat(klon))
372	do ig0=1,klon
373	j=index_j(ig0)
374	klat(ig0)=j
375	enddo
376	endif ! planet_type=titan
377	c------------------------------------------------------------------
378
379	c$OMP END MASTER
380	c$OMP BARRIER
381	ELSE
382	debut = .FALSE.
383	ENDIF
384
385
386	c-----------------------------------------------------------------------
387	c 40. transformation des variables dynamiques en variables physiques:
388	c ---------------------------------------------------------------
389
390	c 41. pressions au sol (en Pascals)
391	c ----------------------------------
392
393	c$OMP MASTER
394	call start_timer(timer_physic)
395	c$OMP END MASTER
396
397	c$OMP MASTER
398	!CDIR ON_ADB(index_i)
399	!CDIR ON_ADB(index_j)
400	do ig0=1,klon
401	i=index_i(ig0)
402	j=index_j(ig0)
403	zpsrf(ig0)=pps(i,j)
404	enddo
405	c$OMP END MASTER
406
407
408	c 42. pression intercouches et fonction d'Exner:
409
410	c -----------------------------------------------------------------
411	c .... zplev definis aux (llm +1) interfaces des couches ....
412	c .... zplay definis aux ( llm ) milieux des couches ....
413	c -----------------------------------------------------------------
414
415	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
416
417	c print *,omp_rank,'klon--->',klon
418	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
419	DO l = 1, llmp1
420	!CDIR ON_ADB(index_i)
421	!CDIR ON_ADB(index_j)
422	do ig0=1,klon
423	i=index_i(ig0)
424	j=index_j(ig0)
425	zplev( ig0,l ) = pp(i,j,l)
426	enddo
427	ENDDO
428	c$OMP END DO NOWAIT
429	if (flag_moyzon) then
430	call AllGather_Field(pp,iip1*jjp1,llmp1)
431	j=index_j(1)
432	tmpvarp1(:,:) = pp(:,j,:)
433	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
434	zplevbar_mpi(1,:) = tmpvarbarp1
435	do ig0=2,klon
436	j=index_j(ig0)
437	if (j.ne.index_j(ig0-1)) then
438	tmpvarp1(:,:) = pp(:,j,:)
439	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
440	zplevbar_mpi(ig0,:) = tmpvarbarp1
441	else
442	zplevbar_mpi(ig0,:) = zplevbar_mpi(ig0-1,:)
443	endif
444	enddo
445	endif
446
447	! ADAPTATION GCM POUR CP(T)
448	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
449	DO l=1,llm
450	!CDIR ON_ADB(index_i)
451	!CDIR ON_ADB(index_j)
452	do ig0=1,klon
453	i=index_i(ig0)
454	j=index_j(ig0)
455	zpk(ig0,l)=ppk(i,j,l)
456	zteta(ig0,l)=pteta(i,j,l)
457	enddo
458	ENDDO
459	c$OMP END DO NOWAIT
460	if (flag_moyzon) then
461	call AllGather_Field(pteta,iip1*jjp1,llm)
462	call AllGather_Field(ppk,iip1*jjp1,llm)
463	j=index_j(1)
464	tmpvar(:,:) = pteta(:,j,:)
465	call moyzon(llm,tmpvar,tmpvarbar)
466	ztetabar_mpi(1,:) = tmpvarbar
467	tmpvar(:,:) = ppk(:,j,:)
468	call moyzon(llm,tmpvar,tmpvarbar)
469	zpkbar_mpi(1,:) = tmpvarbar
470	call tpot2t_p(1,llm,ztetabar_mpi(1,:),ztfibar_mpi(1,:),
471	& zpkbar_mpi(1,:))
472	do ig0=2,klon
473	j=index_j(ig0)
474	if (j.ne.index_j(ig0-1)) then
475	tmpvar(:,:) = pteta(:,j,:)
476	call moyzon(llm,tmpvar,tmpvarbar)
477	ztetabar_mpi(ig0,:) = tmpvarbar
478	tmpvar(:,:) = ppk(:,j,:)
479	call moyzon(llm,tmpvar,tmpvarbar)
480	zpkbar_mpi(ig0,:) = tmpvarbar
481	call tpot2t_p(1,llm,ztetabar_mpi(ig0,:),ztfibar_mpi(ig0,:),
482	& zpkbar_mpi(ig0,:))
483	else
484	zpkbar_mpi(ig0,:) = zpkbar_mpi(ig0-1,:)
485	ztetabar_mpi(ig0,:) = ztetabar_mpi(ig0-1,:)
486	ztfibar_mpi(ig0,:) = ztfibar_mpi(ig0-1,:)
487	endif
488	enddo
489	endif
490
491	c 43. temperature naturelle (en K) et pressions milieux couches .
492	c ---------------------------------------------------------------
493
494	! ADAPTATION GCM POUR CP(T)
495	call tpot2t_p(klon,llm,zteta,ztfi,zpk)
496
497	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
498	DO l=1,llm
499	!CDIR ON_ADB(index_i)
500	!CDIR ON_ADB(index_j)
501	do ig0=1,klon
502	i=index_i(ig0)
503	j=index_j(ig0)
504	pksurcp = ppk(i,j,l) / cpp
505	zplay(ig0,l) = preff * pksurcp ** unskap
506	enddo
507	ENDDO
508	c$OMP END DO NOWAIT
509	if (flag_moyzon) then
510	zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap
511	endif
512
513	c 43.bis traceurs (tous intensifs)
514	c ---------------
515
516	DO iq=1,nqtot
517	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
518	DO l=1,llm
519	!CDIR ON_ADB(index_i)
520	!CDIR ON_ADB(index_j)
521	do ig0=1,klon
522	i=index_i(ig0)
523	j=index_j(ig0)
524	zqfi(ig0,l,iq) = pq(i,j,l,iq)
525	enddo
526	ENDDO
527	c$OMP END DO NOWAIT
528	ENDDO ! of DO iq=1,nqtot
529	if (flag_moyzon) then
530	DO iq=1,nqtot
531	call AllGather_Field(pq(:,:,:,iq),iip1*jjp1,llm)
532	j=index_j(1)
533	tmpvar(:,:) = pq(:,j,:,iq)
534	call moyzon(llm,tmpvar,tmpvarbar)
535	zqfibar_mpi(1,:,iq) = tmpvarbar
536	do ig0=2,klon
537	j=index_j(ig0)
538	if (j.ne.index_j(ig0-1)) then
539	tmpvar(:,:) = pq(:,j,:,iq)
540	call moyzon(llm,tmpvar,tmpvarbar)
541	zqfibar_mpi(ig0,:,iq) = tmpvarbar
542	else
543	zqfibar_mpi(ig0,:,iq) = zqfibar_mpi(ig0-1,:,iq)
544	endif
545	enddo
546	ENDDO ! of DO iq=1,nqtot
547	endif
548
549
550	c Geopotentiel calcule par rapport a la surface locale:
551	c -----------------------------------------------------
552
553	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
554
555	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
556
557	c$OMP BARRIER
558	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
559	DO l=1,llm
560	DO ig=1,klon
561	zphi(ig,l)=zphi(ig,l)-zphis(ig)
562	ENDDO
563	ENDDO
564	c$OMP END DO NOWAIT
565
566	if (flag_moyzon) then
567	call AllGather_Field(pphis,iip1*jjp1,1)
568	call AllGather_Field(pphi,iip1*jjp1,llm)
569	j=index_j(1)
570	tmpvar(:,1) = pphis(:,j)
571	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
572	zphisbar_mpi(1) = tmpvarbar(1)
573	tmpvar(:,:) = pphi(:,j,:)
574	call moyzon(llm,tmpvar,tmpvarbar)
575	zphibar_mpi(1,:) = tmpvarbar-zphisbar_mpi(1)
576	do ig0=2,klon
577	j=index_j(ig0)
578	if (j.ne.index_j(ig0-1)) then
579	tmpvar(:,1) = pphis(:,j)
580	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
581	zphisbar_mpi(ig0) = tmpvarbar(1)
582	tmpvar(:,:) = pphi(:,j,:)
583	call moyzon(llm,tmpvar,tmpvarbar)
584	zphibar_mpi(ig0,:) = tmpvarbar-zphisbar_mpi(ig0)
585	else
586	zphisbar_mpi(ig0) = zphisbar_mpi(ig0-1)
587	zphibar_mpi(ig0,:) = zphibar_mpi(ig0-1,:)
588	endif
589	enddo
590	endif
591
592	c
593	c 45. champ u:
594	c ------------
595
596	kstart=1
597	kend=klon
598
599	if (is_north_pole_dyn) kstart=2
600	if (is_south_pole_dyn) kend=klon-1
601
602	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
603	DO l=1,llm
604	!CDIR ON_ADB(index_i)
605	!CDIR ON_ADB(index_j)
606	!CDIR SPARSE
607	do ig0=kstart,kend
608	i=index_i(ig0)
609	j=index_j(ig0)
610	if (i==1) then
611	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
612	$ + pucov(1,j,l)/cu(1,j) )
613	else
614	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
615	$ + pucov(i,j,l)/cu(i,j) )
616	endif
617	enddo
618	ENDDO
619	c$OMP END DO NOWAIT
620
621	c
622	C Alvaro de la Camara (May 2014)
623	C 46.1 Calcul de la vorticite et passage sur la grille physique
624	C --------------------------------------------------------------
625
626	jjb=jj_begin_dyn-1
627	jje=jj_end_dyn+1
628	if (is_north_pole_dyn) jjb=1
629	if (is_south_pole_dyn) jje=jjm
630
631	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
632
633	DO l=1,llm
634	do i=1,iim
635	do j=jjb,jje
636	zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
637	$ + pucov(i,j+1,l) - pucov(i,j,l))
638	$ / (cu(i,j)+cu(i,j+1))
639	$ / (cv(i+1,j)+cv(i,j)) *4
640	enddo
641	enddo
642	ENDDO
643
644
645	c 46.2champ v:
646	c -----------
647
648	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
649	DO l=1,llm
650	!CDIR ON_ADB(index_i)
651	!CDIR ON_ADB(index_j)
652	DO ig0=kstart,kend
653	i=index_i(ig0)
654	j=index_j(ig0)
655	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
656	$ + pvcov(i,j,l)/cv(i,j) )
657
658	if (j==1 .OR. j==jjp1) then ! AdlC MAY 2014
659	zrfi(ig0,l) = 0 ! AdlC MAY 2014
660	else
661	if(i==1)then
662	zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
663	$ +zrot(1,j-1,l)+zrot(1,j,l)) ! AdlC MAY 2014
664	else
665	zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
666	$ +zrot(i,j-1,l)+zrot(i,j,l)) ! AdlC MAY 2014
667	endif
668	endif
669
670
671	ENDDO
672	ENDDO
673	c$OMP END DO NOWAIT
674
675	c 47. champs de vents aux pole nord
676	c ------------------------------
677	c U = 1 / pi * integrale [ v * cos(long) * d long ]
678	c V = 1 / pi * integrale [ v * sin(long) * d long ]
679
680	if (is_north_pole_dyn) then
681	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
682	DO l=1,llm
683
684	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
685	DO i=2,iim
686	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
687	ENDDO
688
689	DO i=1,iim
690	zcos(i) = COS(rlonv(i))*z1(i)
691	zsin(i) = SIN(rlonv(i))*z1(i)
692	ENDDO
693
694	zufi(1,l) = SSUM(iim,zcos,1)/pi
695	zvfi(1,l) = SSUM(iim,zsin,1)/pi
696	zrfi(1,l) = 0.
697
698	ENDDO
699	c$OMP END DO NOWAIT
700	endif
701
702
703	c 48. champs de vents aux pole sud:
704	c ---------------------------------
705	c U = 1 / pi * integrale [ v * cos(long) * d long ]
706	c V = 1 / pi * integrale [ v * sin(long) * d long ]
707
708	if (is_south_pole_dyn) then
709	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
710	DO l=1,llm
711
712	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
713	DO i=2,iim
714	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
715	ENDDO
716
717	DO i=1,iim
718	zcos(i) = COS(rlonv(i))*z1(i)
719	zsin(i) = SIN(rlonv(i))*z1(i)
720	ENDDO
721
722	zufi(klon,l) = SSUM(iim,zcos,1)/pi
723	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
724	zrfi(klon,l) = 0.
725	ENDDO
726	c$OMP END DO NOWAIT
727	endif
728
729	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
730	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
731
732	c-----------------------------------------------------------------------
733	c Appel de la physique:
734	c ---------------------
735
736	! Appel de la physique: pose probleme quand on tourne
737	! SANS physique, car physiq.F est dans le repertoire phy[]...
738	! Il faut une cle CPP_PHYS
739
740	! Le fait que les arguments de physiq soient differents selon les planetes
741	! ne pose pas de probleme a priori.
742
743
744	c$OMP BARRIER
745	if (first_omp) then
746	klon=klon_omp
747
748	allocate(zplev_omp(klon,llm+1))
749	allocate(zplay_omp(klon,llm))
750	allocate(zpk_omp(klon,llm))
751	allocate(zphi_omp(klon,llm))
752	allocate(zphis_omp(klon))
753	allocate(presnivs_omp(llm))
754	allocate(zufi_omp(klon,llm))
755	allocate(zvfi_omp(klon,llm))
756	allocate(zrfi_omp(klon,llm)) ! LG Ari 2014
757	allocate(ztfi_omp(klon,llm))
758	allocate(zqfi_omp(klon,llm,nqtot))
759	allocate(zdufi_omp(klon,llm))
760	allocate(zdvfi_omp(klon,llm))
761	allocate(zdtfi_omp(klon,llm))
762	allocate(zdqfi_omp(klon,llm,nqtot))
763	allocate(zdufic_omp(klon,llm))
764	allocate(zdvfic_omp(klon,llm))
765	allocate(zdtfic_omp(klon,llm))
766	allocate(zdqfic_omp(klon,llm,nqtot))
767	allocate(zdpsrf_omp(klon))
768	allocate(flxwfi_omp(klon,llm))
769
770	if (flag_moyzon) call moyzon_init_omp(klon,llm,nqtot)
771
772	first_omp=.false.
773	endif
774
775
776	klon=klon_omp
777	offset=klon_omp_begin-1
778
779	do l=1,llm+1
780	do i=1,klon
781	zplev_omp(i,l)=zplev(offset+i,l)
782	enddo
783	enddo
784
785	do l=1,llm
786	do i=1,klon
787	zplay_omp(i,l)=zplay(offset+i,l)
788	enddo
789	enddo
790
791	do l=1,llm
792	do i=1,klon
793	zpk_omp(i,l)=zpk(offset+i,l)
794	enddo
795	enddo
796
797	do l=1,llm
798	do i=1,klon
799	zphi_omp(i,l)=zphi(offset+i,l)
800	enddo
801	enddo
802
803	do i=1,klon
804	zphis_omp(i)=zphis(offset+i)
805	enddo
806
807
808	do l=1,llm
809	presnivs_omp(l)=presnivs(l)
810	enddo
811
812	do l=1,llm
813	do i=1,klon
814	zufi_omp(i,l)=zufi(offset+i,l)
815	enddo
816	enddo
817
818	do l=1,llm
819	do i=1,klon
820	zvfi_omp(i,l)=zvfi(offset+i,l)
821	enddo
822	enddo
823
824	do l=1,llm
825	do i=1,klon
826	zrfi_omp(i,l)=zrfi(offset+i,l)
827	enddo
828	enddo
829
830
831	do l=1,llm
832	do i=1,klon
833	ztfi_omp(i,l)=ztfi(offset+i,l)
834	enddo
835	enddo
836
837	do iq=1,nqtot
838	do l=1,llm
839	do i=1,klon
840	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
841	enddo
842	enddo
843	enddo
844
845	do l=1,llm
846	do i=1,klon
847	zdufi_omp(i,l)=zdufi(offset+i,l)
848	enddo
849	enddo
850
851	do l=1,llm
852	do i=1,klon
853	zdvfi_omp(i,l)=zdvfi(offset+i,l)
854	enddo
855	enddo
856
857	do l=1,llm
858	do i=1,klon
859	zdtfi_omp(i,l)=zdtfi(offset+i,l)
860	enddo
861	enddo
862
863	do iq=1,nqtot
864	do l=1,llm
865	do i=1,klon
866	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
867	enddo
868	enddo
869	enddo
870
871	do i=1,klon
872	zdpsrf_omp(i)=zdpsrf(offset+i)
873	enddo
874
875	do l=1,llm
876	do i=1,klon
877	flxwfi_omp(i,l)=flxwfi(offset+i,l)
878	enddo
879	enddo
880
881	if (flag_moyzon) then
882	do l=1,llm+1
883	do i=1,klon
884	zplevbar(i,l)=zplevbar_mpi(offset+i,l)
885	enddo
886	enddo
887
888	do l=1,llm
889	do i=1,klon
890	zplaybar(i,l)=zplaybar_mpi(offset+i,l)
891	enddo
892	enddo
893
894	do l=1,llm
895	do i=1,klon
896	zphibar(i,l)=zphibar_mpi(offset+i,l)
897	enddo
898	enddo
899
900	do i=1,klon
901	zphisbar(i)=zphisbar_mpi(offset+i)
902	enddo
903
904	do l=1,llm
905	do i=1,klon
906	ztfibar(i,l)=ztfibar_mpi(offset+i,l)
907	enddo
908	enddo
909
910	do iq=1,nqtot
911	do l=1,llm
912	do i=1,klon
913	zqfibar(i,l,iq)=zqfibar_mpi(offset+i,l,iq)
914	enddo
915	enddo
916	enddo
917	endif
918
919
920	c$OMP BARRIER
921
922	!$OMP MASTER
923	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
924	!$OMP END MASTER
925	zdt_split=dtphys/nsplit_phys
926	zdufic_omp(:,:)=0.
927	zdvfic_omp(:,:)=0.
928	zdtfic_omp(:,:)=0.
929	zdqfic_omp(:,:,:)=0.
930
931	#ifdef CPP_PHYS
932	do isplit=1,nsplit_phys
933
934	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
935	debut_split=debut.and.isplit==1
936	lafin_split=lafin.and.isplit==nsplit_phys
937
938	CALL call_physiq(klon,llm,nqtot,tname,
939	& debut_split,lafin_split,
940	& jD_cur,jH_cur_split,zdt_split,
941	& zplev_omp,zplay_omp,
942	& zpk_omp,zphi_omp,zphis_omp,
943	& presnivs_omp,
944	& zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
945	& flxwfi_omp,pducov,
946	& zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
947	& zdpsrf_omp,tracerdyn)
948
949	! if (planet_type=="earth") then
950	! CALL physiq (klon,
951	! . llm,
952	! . debut_split,
953	! . lafin_split,
954	! . jD_cur,
955	! . jH_cur_split,
956	! . zdt_split,
957	! . zplev_omp,
958	! . zplay_omp,
959	! . zphi_omp,
960	! . zphis_omp,
961	! . presnivs_omp,
962	! . zufi_omp,
963	! . zvfi_omp,
964	! . ztfi_omp,
965	! . zqfi_omp,
966	! . flxwfi_omp,
967	! . zdufi_omp,
968	! . zdvfi_omp,
969	! . zdtfi_omp,
970	! . zdqfi_omp,
971	! . zdpsrf_omp,
972	! . pducov)
973	!
974	! else if ( planet_type=="generic" ) then
975	!
976	! CALL physiq (klon, !! ngrid
977	! . llm, !! nlayer
978	! . nqtot, !! nq
979	! . tname, !! tracer names from dynamical core (given in infotrac)
980	! . debut_split, !! firstcall
981	! . lafin_split, !! lastcall
982	! . jD_cur, !! pday. see leapfrog_p
983	! . jH_cur_split, !! ptime "fraction of day"
984	! . zdt_split, !! ptimestep
985	! . zplev_omp, !! pplev
986	! . zplay_omp, !! pplay
987	! . zphi_omp, !! pphi
988	! . zufi_omp, !! pu
989	! . zvfi_omp, !! pv
990	! . ztfi_omp, !! pt
991	! . zqfi_omp, !! pq
992	! . flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
993	! . zdufi_omp, !! pdu
994	! . zdvfi_omp, !! pdv
995	! . zdtfi_omp, !! pdt
996	! . zdqfi_omp, !! pdq
997	! . zdpsrf_omp, !! pdpsrf
998	! . tracerdyn) !! tracerdyn <-- utilite ???
999	!
1000	! else if ( planet_type=="mars" ) then
1001	!
1002	! CALL physiq (klon, ! ngrid
1003	! . llm, ! nlayer
1004	! . nqtot, ! nq
1005	! . debut_split, ! firstcall
1006	! . lafin_split, ! lastcall
1007	! . jD_cur, ! pday
1008	! . jH_cur_split, ! ptime
1009	! . zdt_split, ! ptimestep
1010	! . zplev_omp, ! pplev
1011	! . zplay_omp, ! pplay
1012	! . zphi_omp, ! pphi
1013	! . zufi_omp, ! pu
1014	! . zvfi_omp, ! pv
1015	! . ztfi_omp, ! pt
1016	! . zqfi_omp, ! pq
1017	! . flxwfi_omp, ! pw
1018	! . zdufi_omp, ! pdu
1019	! . zdvfi_omp, ! pdv
1020	! . zdtfi_omp, ! pdt
1021	! . zdqfi_omp, ! pdq
1022	! . zdpsrf_omp, ! pdpsrf
1023	! . tracerdyn) ! tracerdyn (somewhat obsolete)
1024	!
1025	! else if ((planet_type=="titan").or.(planet_type=="venus")) then
1026	!
1027	! CALL physiq (klon,
1028	! . llm,
1029	! . nqtot,
1030	! . debut_split,
1031	! . lafin_split,
1032	! . jD_cur,
1033	! . jH_cur_split,
1034	! . zdt_split,
1035	! . zplev_omp,
1036	! . zplay_omp,
1037	! . zpk_omp,
1038	! . zphi_omp,
1039	! . zphis_omp,
1040	! . presnivs_omp,
1041	! . zufi_omp,
1042	! . zvfi_omp,
1043	! . ztfi_omp,
1044	! . zqfi_omp,
1045	! . flxwfi_omp,
1046	! . zdufi_omp,
1047	! . zdvfi_omp,
1048	! . zdtfi_omp,
1049	! . zdqfi_omp,
1050	! . zdpsrf_omp)
1051	!
1052	! else ! unknown "planet_type"
1053	!
1054	! write(lunout,*) "calfis_p: error, unknown planet_type: ",
1055	! & trim(planet_type)
1056	! stop
1057	!
1058	! endif ! planet_type
1059	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
1060	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
1061	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
1062	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
1063
1064	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
1065	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
1066	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
1067	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
1068
1069	enddo
1070
1071	#endif
1072	! of #ifdef CPP_PHYS
1073
1074	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
1075	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
1076	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
1077	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
1078
1079	c$OMP BARRIER
1080
1081	do l=1,llm+1
1082	do i=1,klon
1083	zplev(offset+i,l)=zplev_omp(i,l)
1084	enddo
1085	enddo
1086
1087	do l=1,llm
1088	do i=1,klon
1089	zplay(offset+i,l)=zplay_omp(i,l)
1090	enddo
1091	enddo
1092
1093	do l=1,llm
1094	do i=1,klon
1095	zphi(offset+i,l)=zphi_omp(i,l)
1096	enddo
1097	enddo
1098
1099
1100	do i=1,klon
1101	zphis(offset+i)=zphis_omp(i)
1102	enddo
1103
1104
1105	do l=1,llm
1106	presnivs(l)=presnivs_omp(l)
1107	enddo
1108
1109	do l=1,llm
1110	do i=1,klon
1111	zufi(offset+i,l)=zufi_omp(i,l)
1112	enddo
1113	enddo
1114
1115	do l=1,llm
1116	do i=1,klon
1117	zvfi(offset+i,l)=zvfi_omp(i,l)
1118	enddo
1119	enddo
1120
1121	do l=1,llm
1122	do i=1,klon
1123	ztfi(offset+i,l)=ztfi_omp(i,l)
1124	enddo
1125	enddo
1126
1127	do iq=1,nqtot
1128	do l=1,llm
1129	do i=1,klon
1130	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
1131	enddo
1132	enddo
1133	enddo
1134
1135	do l=1,llm
1136	do i=1,klon
1137	zdufi(offset+i,l)=zdufi_omp(i,l)
1138	enddo
1139	enddo
1140
1141	do l=1,llm
1142	do i=1,klon
1143	zdvfi(offset+i,l)=zdvfi_omp(i,l)
1144	enddo
1145	enddo
1146
1147	do l=1,llm
1148	do i=1,klon
1149	zdtfi(offset+i,l)=zdtfi_omp(i,l)
1150	enddo
1151	enddo
1152
1153	do iq=1,nqtot
1154	do l=1,llm
1155	do i=1,klon
1156	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
1157	enddo
1158	enddo
1159	enddo
1160
1161	do i=1,klon
1162	zdpsrf(offset+i)=zdpsrf_omp(i)
1163	enddo
1164
1165
1166	klon=klon_mpi
1167	500 CONTINUE
1168	c$OMP BARRIER
1169
1170	c$OMP MASTER
1171	call stop_timer(timer_physic)
1172	c$OMP END MASTER
1173
1174	IF (using_mpi) THEN
1175
1176	if (MPI_rank>0) then
1177
1178	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1179	DO l=1,llm
1180	du_send(1:iim,l)=zdufi(1:iim,l)
1181	dv_send(1:iim,l)=zdvfi(1:iim,l)
1182	ENDDO
1183	c$OMP END DO NOWAIT
1184
1185	c$OMP BARRIER
1186	#ifdef CPP_MPI
1187	c$OMP MASTER
1188	!$OMP CRITICAL (MPI)
1189	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
1190	& COMM_LMDZ,Req(1),ierr)
1191	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
1192	& COMM_LMDZ,Req(2),ierr)
1193	!$OMP END CRITICAL (MPI)
1194	c$OMP END MASTER
1195	#endif
1196	c$OMP BARRIER
1197
1198	endif
1199
1200	if (MPI_rank<MPI_Size-1) then
1201	c$OMP BARRIER
1202	#ifdef CPP_MPI
1203	c$OMP MASTER
1204	!$OMP CRITICAL (MPI)
1205	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
1206	& COMM_LMDZ,Req(3),ierr)
1207	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
1208	& COMM_LMDZ,Req(4),ierr)
1209	!$OMP END CRITICAL (MPI)
1210	c$OMP END MASTER
1211	#endif
1212	endif
1213
1214	c$OMP BARRIER
1215
1216
1217	#ifdef CPP_MPI
1218	c$OMP MASTER
1219	!$OMP CRITICAL (MPI)
1220	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
1221	call MPI_WAITALL(4,Req(1),Status,ierr)
1222	else if (MPI_rank>0) then
1223	call MPI_WAITALL(2,Req(1),Status,ierr)
1224	else if (MPI_rank <MPI_Size-1) then
1225	call MPI_WAITALL(2,Req(3),Status,ierr)
1226	endif
1227	!$OMP END CRITICAL (MPI)
1228	c$OMP END MASTER
1229	#endif
1230
1231	c$OMP BARRIER
1232
1233	ENDIF ! using_mpi
1234
1235
1236	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1237	DO l=1,llm
1238
1239	zdufi2(1:klon,l)=zdufi(1:klon,l)
1240	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
1241
1242	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
1243	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
1244
1245	pdhfi(:,jj_begin,l)=0
1246	pdqfi(:,jj_begin,l,:)=0
1247	pdufi(:,jj_begin,l)=0
1248	pdvfi(:,jj_begin,l)=0
1249
1250	if (.not. is_south_pole_dyn) then
1251	pdhfi(:,jj_end,l)=0
1252	pdqfi(:,jj_end,l,:)=0
1253	pdufi(:,jj_end,l)=0
1254	pdvfi(:,jj_end,l)=0
1255	endif
1256
1257	ENDDO
1258	c$OMP END DO NOWAIT
1259
1260	c$OMP MASTER
1261	pdpsfi(:,jj_begin)=0
1262	if (.not. is_south_pole_dyn) then
1263	pdpsfi(:,jj_end)=0
1264	endif
1265	c$OMP END MASTER
1266	c-----------------------------------------------------------------------
1267	c transformation des tendances physiques en tendances dynamiques:
1268	c ---------------------------------------------------------------
1269
1270	c tendance sur la pression :
1271	c -----------------------------------
1272	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
1273	c
1274	c 62. enthalpie potentielle
1275	c ---------------------
1276
1277	kstart=1
1278	kend=klon
1279
1280	if (is_north_pole_dyn) kstart=2
1281	if (is_south_pole_dyn) kend=klon-1
1282
1283	! ADAPTATION GCM POUR CP(T)
1284	call t2tpot_p(klon,llm,ztfi,zteta,zpk)
1285
1286
1287	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1288	DO l=1,llm
1289	!CDIR ON_ADB(index_i)
1290	!CDIR ON_ADB(index_j)
1291	!cdir NODEP
1292	do ig0=kstart,kend
1293	i=index_i(ig0)
1294	j=index_j(ig0)
1295	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1296	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1297	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1298	if (i==1) then
1299	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1300	endif
1301	enddo
1302
1303	if (is_north_pole_dyn) then
1304	DO i=1,iip1
1305	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
1306	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
1307	enddo
1308	endif
1309
1310	if (is_south_pole_dyn) then
1311	DO i=1,iip1
1312	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1313	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
1314	ENDDO
1315	endif
1316	ENDDO
1317	c$OMP END DO NOWAIT
1318
1319	c 62. humidite specifique
1320	c ---------------------
1321	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1322	! DO iq=1,nqtot
1323	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1324	! DO l=1,llm
1325	!!!cdir NODEP
1326	! do ig0=kstart,kend
1327	! i=index_i(ig0)
1328	! j=index_j(ig0)
1329	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1330	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1331	! enddo
1332	!
1333	! if (is_north_pole_dyn) then
1334	! do i=1,iip1
1335	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1336	! enddo
1337	! endif
1338	!
1339	! if (is_south_pole_dyn) then
1340	! do i=1,iip1
1341	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1342	! enddo
1343	! endif
1344	! ENDDO
1345	!c$OMP END DO NOWAIT
1346	! ENDDO
1347
1348	c 63. traceurs (tous en intensifs)
1349	c ------------
1350	C initialisation des tendances
1351
1352	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1353	DO l=1,llm
1354	pdqfi(:,:,l,:)=0.
1355	ENDDO
1356	c$OMP END DO NOWAIT
1357
1358	C
1359	!cdir NODEP
1360	DO iq=1,nqtot
1361	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1362	DO l=1,llm
1363	!CDIR ON_ADB(index_i)
1364	!CDIR ON_ADB(index_j)
1365	!cdir NODEP
1366	DO ig0=kstart,kend
1367	i=index_i(ig0)
1368	j=index_j(ig0)
1369	pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1370	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1371	ENDDO
1372
1373	IF (is_north_pole_dyn) then
1374	DO i=1,iip1
1375	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1376	ENDDO
1377	ENDIF
1378
1379	IF (is_south_pole_dyn) then
1380	DO i=1,iip1
1381	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1382	ENDDO
1383	ENDIF
1384
1385	ENDDO
1386	c$OMP END DO NOWAIT
1387	ENDDO
1388
1389	c 65. champ u:
1390	c ------------
1391	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1392	DO l=1,llm
1393	!CDIR ON_ADB(index_i)
1394	!CDIR ON_ADB(index_j)
1395	!cdir NODEP
1396	do ig0=kstart,kend
1397	i=index_i(ig0)
1398	j=index_j(ig0)
1399
1400	if (i/=iim) then
1401	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1402	endif
1403
1404	if (i==1) then
1405	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1406	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1407	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1408	endif
1409
1410	enddo
1411
1412	if (is_north_pole_dyn) then
1413	DO i=1,iip1
1414	pdufi(i,1,l) = 0.
1415	ENDDO
1416	endif
1417
1418	if (is_south_pole_dyn) then
1419	DO i=1,iip1
1420	pdufi(i,jjp1,l) = 0.
1421	ENDDO
1422	endif
1423
1424	ENDDO
1425	c$OMP END DO NOWAIT
1426
1427	c 67. champ v:
1428	c ------------
1429
1430	kstart=1
1431	kend=klon
1432
1433	if (is_north_pole_dyn) kstart=2
1434	if (is_south_pole_dyn) kend=klon-1-iim
1435
1436	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1437	DO l=1,llm
1438	!CDIR ON_ADB(index_i)
1439	!CDIR ON_ADB(index_j)
1440	!cdir NODEP
1441	do ig0=kstart,kend
1442	i=index_i(ig0)
1443	j=index_j(ig0)
1444	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1445	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1446	$ zdvfi2(ig0+iim,l))
1447	$ *cv(i,j)
1448	enddo
1449
1450	ENDDO
1451	c$OMP END DO NOWAIT
1452
1453
1454	c 68. champ v pres des poles:
1455	c ---------------------------
1456	c v = U * cos(long) + V * SIN(long)
1457
1458	if (is_north_pole_dyn) then
1459
1460	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1461	DO l=1,llm
1462
1463	DO i=1,iim
1464	pdvfi(i,1,l)=
1465	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1466
1467	pdvfi(i,1,l)=
1468	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1469	ENDDO
1470
1471	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1472
1473	ENDDO
1474	c$OMP END DO NOWAIT
1475
1476	endif
1477
1478	if (is_south_pole_dyn) then
1479
1480	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1481	DO l=1,llm
1482
1483	DO i=1,iim
1484	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1485	$ +zdvfi(klon,l)*SIN(rlonv(i))
1486
1487	pdvfi(i,jjm,l)=
1488	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1489	ENDDO
1490
1491	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1492
1493	ENDDO
1494	c$OMP END DO NOWAIT
1495
1496	endif
1497	c-----------------------------------------------------------------------
1498
1499	700 CONTINUE
1500
1501	firstcal = .FALSE.
1502
1503	#else
1504	write(lunout,*)
1505	& "calfis_p: for now can only work with parallel physics"
1506	stop
1507	#endif
1508	! of #ifdef CPP_PHYS
1509	#endif
1510	! of #ifdef CPP_PARA
1511	END

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

Download in other formats:

Original Format