Context Navigation

calfis_p.F @ 1543

Last change on this file since 1543 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

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

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

Download in other formats:

Original Format