source: trunk/LMDZ.COMMON/libf/dynphy_lonlat/calfis_p.F @ 3553

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

All models:
Clean up the dynamics/physics interface to converge with LMDZ5;
get rid of tracerdyn parameter (which was supposed to send information
about wether tracers should be advected or not in the Mars and Generic
models).
EM

File size: 42.4 KB
Line 
1!
2! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fairhead $
3!
4C
5C
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
58c=======================================================================
59c
60c   1. rearrangement des tableaux et transformation
61c      variables dynamiques  >  variables physiques
62c   2. calcul des termes physiques
63c   3. retransformation des tendances physiques en tendances dynamiques
64c
65c   remarques:
66c   ----------
67c
68c    - les vents sont donnes dans la physique par leurs composantes
69c      naturelles.
70c    - la variable thermodynamique de la physique est une variable
71c      intensive :   T
72c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
73c    - les deux seules variables dependant de la geometrie necessaires
74c      pour la physique sont la latitude pour le rayonnement et
75c      l'aire de la maille quand on veut integrer une grandeur
76c      horizontalement.
77c    - les points de la physique sont les points scalaires de la
78c      la dynamique; numerotation:
79c          1 pour le pole nord
80c          (jjm-1)*iim pour l'interieur du domaine
81c          ngridmx pour le pole sud
82c      ---> ngridmx=2+(jjm-1)*iim
83c
84c     Input :
85c     -------
86c       ecritphy        frequence d'ecriture (en jours)de histphy
87c       pucov           covariant zonal velocity
88c       pvcov           covariant meridional velocity
89c       pteta           potential temperature
90c       pps             surface pressure
91c       pmasse          masse d'air dans chaque maille
92c       pts             surface temperature  (K)
93c       callrad         clef d'appel au rayonnement
94c
95c    Output :
96c    --------
97c        pdufi          tendency for the natural zonal velocity (ms-1)
98c        pdvfi          tendency for the natural meridional velocity
99c        pdhfi          tendency for the potential temperature (K/s)
100c        pdtsfi         tendency for the surface temperature
101c
102c        pdtrad         radiative tendencies  \  both input
103c        pfluxrad       radiative fluxes      /  and output
104c
105c=======================================================================
106c
107c-----------------------------------------------------------------------
108c
109c    0.  Declarations :
110c    ------------------
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
123c    Arguments :
124c    -----------
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
158c    Local variables :
159c    -----------------
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(:,:,:)
178c
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
224c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp,
225c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
226c$OMP+                 zrfi_omp,zqfi_omp,zdufi_omp,zdvfi_omp,
227c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
228c$OMP+                 zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)       
229
230      LOGICAL,SAVE :: first_omp=.true.
231c$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.
241c$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
265c-----------------------------------------------------------------------
266
267c    1. Initialisations :
268c    --------------------
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
286c$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
315c------------------------------------------------------------------
316c 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)
358c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
359        zlaymoy(1:llm) = g*rad*rad/(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)=(zz1*zlaymoy(l-1)+zz2*zlaymoy(l))/(zz1+zz2)
365        ENDDO
366        zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
367c-------------------
368c + 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
375c------------------------------------------------------------------
376
377c$OMP END MASTER
378c$OMP BARRIER     
379      ELSE
380          debut = .FALSE.
381      ENDIF
382
383
384c-----------------------------------------------------------------------
385c   40. transformation des variables dynamiques en variables physiques:
386c   ---------------------------------------------------------------
387
388c   41. pressions au sol (en Pascals)
389c   ----------------------------------
390
391c$OMP MASTER
392      call start_timer(timer_physic)
393c$OMP END MASTER
394
395c$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
403c$OMP END MASTER
404
405
406c   42. pression intercouches et fonction d'Exner:
407
408c   -----------------------------------------------------------------
409c     .... zplev  definis aux (llm +1) interfaces des couches  ....
410c     .... zplay  definis aux (  llm )    milieux des couches  ....
411c   -----------------------------------------------------------------
412
413c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
414
415c      print *,omp_rank,'klon--->',klon
416c$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
426c$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)
446c$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
457c$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
489c   43. temperature naturelle (en K) et pressions milieux couches .
490c   ---------------------------------------------------------------
491
492! ADAPTATION GCM POUR CP(T)
493         call tpot2t_p(klon,llm,zteta,ztfi,zpk)
494
495c$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
506c$OMP END DO NOWAIT
507      if (flag_moyzon) then
508        zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap
509      endif
510
511c   43.bis traceurs (tous intensifs)
512c   ---------------
513
514      DO iq=1,nqtot
515c$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
525c$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
548c   Geopotentiel calcule par rapport a la surface locale:
549c   -----------------------------------------------------
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
555c$OMP BARRIER
556c$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
562c$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
590c
591c   45. champ u:
592c   ------------
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     
600c$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
617c$OMP END DO NOWAIT
618
619c
620C  Alvaro de la Camara (May 2014)
621C  46.1 Calcul de la vorticite et passage sur la grille physique
622C  --------------------------------------------------------------
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
629c$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
643c   46.2champ v:
644c   -----------
645
646c$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
671c$OMP END DO NOWAIT
672
673c   47. champs de vents aux pole nord   
674c   ------------------------------
675c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
676c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
677
678      if (is_north_pole_dyn) then
679c$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
697c$OMP END DO NOWAIT     
698      endif
699
700
701c   48. champs de vents aux pole sud:
702c   ---------------------------------
703c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
704c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
705
706      if (is_south_pole_dyn) then
707c$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
724c$OMP END DO NOWAIT       
725      endif
726
727c 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
730c-----------------------------------------------------------------------
731c   Appel de la physique:
732c   ---------------------
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
742c$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     
918c$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
1077c$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
1165500   CONTINUE
1166c$OMP BARRIER
1167
1168c$OMP MASTER
1169      call stop_timer(timer_physic)
1170c$OMP END MASTER
1171
1172      IF (using_mpi) THEN
1173           
1174      if (MPI_rank>0) then
1175
1176c$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
1181c$OMP END DO NOWAIT       
1182
1183c$OMP BARRIER
1184#ifdef CPP_MPI
1185c$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)
1192c$OMP END MASTER
1193#endif
1194c$OMP BARRIER
1195     
1196      endif
1197   
1198      if (MPI_rank<MPI_Size-1) then
1199c$OMP BARRIER
1200#ifdef CPP_MPI
1201c$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)
1208c$OMP END MASTER
1209#endif
1210      endif
1211
1212c$OMP BARRIER
1213
1214
1215#ifdef CPP_MPI
1216c$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)
1226c$OMP END MASTER
1227#endif
1228
1229c$OMP BARRIER     
1230
1231      ENDIF ! using_mpi
1232     
1233     
1234c$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
1256c$OMP END DO NOWAIT
1257
1258c$OMP MASTER
1259       pdpsfi(:,jj_begin)=0   
1260       if (.not. is_south_pole_dyn) then
1261         pdpsfi(:,jj_end)=0
1262       endif
1263c$OMP END MASTER
1264c-----------------------------------------------------------------------
1265c   transformation des tendances physiques en tendances dynamiques:
1266c   ---------------------------------------------------------------
1267
1268c  tendance sur la pression :
1269c  -----------------------------------
1270      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
1271c
1272c   62. enthalpie potentielle
1273c   ---------------------
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
1285c$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
1315c$OMP END DO NOWAIT
1316     
1317c   62. humidite specifique
1318c   ---------------------
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
1346c   63. traceurs (tous en intensifs)
1347c   ------------
1348C     initialisation des tendances
1349
1350c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1351      DO l=1,llm
1352        pdqfi(:,:,l,:)=0.
1353      ENDDO
1354c$OMP END DO NOWAIT     
1355
1356C
1357!cdir NODEP
1358      DO iq=1,nqtot
1359c$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
1384c$OMP END DO NOWAIT     
1385      ENDDO
1386     
1387c   65. champ u:
1388c   ------------
1389c$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
1423c$OMP END DO NOWAIT
1424
1425c   67. champ v:
1426c   ------------
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     
1434c$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
1449c$OMP END DO NOWAIT
1450
1451
1452c   68. champ v pres des poles:
1453c   ---------------------------
1454c      v = U * cos(long) + V * SIN(long)
1455
1456      if (is_north_pole_dyn) then
1457
1458c$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
1472c$OMP END DO NOWAIT
1473
1474      endif   
1475     
1476      if (is_south_pole_dyn) then
1477
1478c$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
1492c$OMP END DO NOWAIT
1493     
1494      endif
1495c-----------------------------------------------------------------------
1496
1497700   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.