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

Last change on this file since 3891 was 3891, checked in by emillour, 4 months ago

Common Dynamics-Physics interface:
Fix a buggy (wrt OpenMP) initialization I introduced in r3825.
EM
File size: 42.7 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! Initialize zphi,zphis to 0 , not really necessary but otherwise
551! debug mode sometimes spots illicit operations on these arrays
552! (in sections actually not used later on, so not really a bug)
553c$OMP MASTER
554      zphi(:,:)=0
555      zphis(:)=0
556c$OMP END MASTER
557c$OMP BARRIER
558     
559      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
560
561      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
562
563c$OMP BARRIER
564c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
565      DO l=1,llm
566         DO ig=1,klon
567           zphi(ig,l)=zphi(ig,l)-zphis(ig)
568         ENDDO
569      ENDDO
570c$OMP END DO NOWAIT
571     
572      if (flag_moyzon) then
573        call AllGather_Field(pphis,iip1*jjp1,1)
574        call AllGather_Field(pphi,iip1*jjp1,llm)
575        j=index_j(1)
576        tmpvar(:,1) = pphis(:,j)
577        call moyzon(1,tmpvar(:,1),tmpvarbar(1))
578        zphisbar_mpi(1) = tmpvarbar(1)
579        tmpvar(:,:) = pphi(:,j,:)
580        call moyzon(llm,tmpvar,tmpvarbar)
581        zphibar_mpi(1,:) = tmpvarbar-zphisbar_mpi(1)
582        do ig0=2,klon
583          j=index_j(ig0)
584          if (j.ne.index_j(ig0-1)) then
585            tmpvar(:,1) = pphis(:,j)
586            call moyzon(1,tmpvar(:,1),tmpvarbar(1))
587            zphisbar_mpi(ig0) = tmpvarbar(1)
588            tmpvar(:,:) = pphi(:,j,:)
589            call moyzon(llm,tmpvar,tmpvarbar)
590            zphibar_mpi(ig0,:) = tmpvarbar-zphisbar_mpi(ig0)
591          else
592            zphisbar_mpi(ig0)  = zphisbar_mpi(ig0-1)
593            zphibar_mpi(ig0,:) = zphibar_mpi(ig0-1,:)
594          endif
595        enddo
596      endif
597
598c
599c   45. champ u:
600c   ------------
601
602      kstart=1
603      kend=klon
604     
605      if (is_north_pole_dyn) kstart=2
606      if (is_south_pole_dyn) kend=klon-1
607     
608c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
609      DO l=1,llm
610!CDIR ON_ADB(index_i)
611!CDIR ON_ADB(index_j)
612!CDIR SPARSE
613        do ig0=kstart,kend
614          i=index_i(ig0)
615          j=index_j(ig0)
616          if (i==1) then
617            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
618     $                         + pucov(1,j,l)/cu(1,j) )
619          else
620            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j)
621     $                       + pucov(i,j,l)/cu(i,j) )
622          endif
623        enddo
624      ENDDO
625c$OMP END DO NOWAIT
626
627c
628C  Alvaro de la Camara (May 2014)
629C  46.1 Calcul de la vorticite et passage sur la grille physique
630C  --------------------------------------------------------------
631
632      jjb=jj_begin_dyn-1
633      jje=jj_end_dyn+1
634      if (is_north_pole_dyn) jjb=1
635      if (is_south_pole_dyn) jje=jjm
636
637c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
638
639      DO l=1,llm
640        do i=1,iim
641          do j=jjb,jje
642            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
643     $                   + pucov(i,j+1,l) - pucov(i,j,l))
644     $                   / (cu(i,j)+cu(i,j+1))
645     $                   / (cv(i+1,j)+cv(i,j)) *4
646          enddo
647        enddo
648      ENDDO
649
650
651c   46.2champ v:
652c   -----------
653
654c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
655      DO l=1,llm
656!CDIR ON_ADB(index_i)
657!CDIR ON_ADB(index_j)
658        DO ig0=kstart,kend
659          i=index_i(ig0)
660          j=index_j(ig0)
661          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1)
662     $                       + pvcov(i,j,l)/cv(i,j) )
663   
664          if (j==1 .OR. j==jjp1) then !  AdlC MAY 2014
665            zrfi(ig0,l) = 0 !  AdlC MAY 2014
666          else
667            if(i==1)then
668            zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
669     $                   +zrot(1,j-1,l)+zrot(1,j,l))   !  AdlC MAY 2014
670            else
671            zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
672     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
673            endif
674          endif
675
676
677         ENDDO
678      ENDDO
679c$OMP END DO NOWAIT
680
681c   47. champs de vents aux pole nord   
682c   ------------------------------
683c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
684c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
685
686      if (is_north_pole_dyn) then
687c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
688        DO l=1,llm
689
690           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
691           DO i=2,iim
692              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
693           ENDDO
694 
695           DO i=1,iim
696              zcos(i)   = COS(rlonv(i))*z1(i)
697              zsin(i)   = SIN(rlonv(i))*z1(i)
698           ENDDO
699 
700           zufi(1,l)  = SSUM(iim,zcos,1)/pi
701           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
702           zrfi(1,l)  = 0.
703 
704        ENDDO
705c$OMP END DO NOWAIT     
706      endif
707
708
709c   48. champs de vents aux pole sud:
710c   ---------------------------------
711c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
712c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
713
714      if (is_south_pole_dyn) then
715c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
716        DO l=1,llm
717 
718         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
719           DO i=2,iim
720             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
721           ENDDO
722 
723           DO i=1,iim
724              zcos(i)    = COS(rlonv(i))*z1(i)
725              zsin(i)    = SIN(rlonv(i))*z1(i)
726           ENDDO
727 
728           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
729           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
730           zrfi(klon,l)  = 0.
731        ENDDO
732c$OMP END DO NOWAIT       
733      endif
734
735c On change de grille, dynamique vers physiq, pour le flux de masse verticale
736      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
737
738c-----------------------------------------------------------------------
739c   Appel de la physique:
740c   ---------------------
741
742! Appel de la physique: pose probleme quand on tourne
743! SANS physique, car physiq.F est dans le repertoire phy[]...
744! Il faut une cle CPP_PHYS
745
746! Le fait que les arguments de physiq soient differents selon les planetes
747! ne pose pas de probleme a priori.
748
749
750c$OMP BARRIER
751      if (first_omp) then
752        klon=klon_omp
753
754        allocate(zplev_omp(klon,llm+1))
755        allocate(zplay_omp(klon,llm))
756        allocate(zpk_omp(klon,llm))
757        allocate(zphi_omp(klon,llm))
758        allocate(zphis_omp(klon))
759        allocate(presnivs_omp(llm))
760        allocate(zufi_omp(klon,llm))
761        allocate(zvfi_omp(klon,llm))
762        allocate(zrfi_omp(klon,llm))  ! LG Ari 2014
763        allocate(ztfi_omp(klon,llm))
764        allocate(zqfi_omp(klon,llm,nqtot))
765        allocate(zdufi_omp(klon,llm))
766        allocate(zdvfi_omp(klon,llm))
767        allocate(zdtfi_omp(klon,llm))
768        allocate(zdqfi_omp(klon,llm,nqtot))
769        allocate(zdufic_omp(klon,llm))
770        allocate(zdvfic_omp(klon,llm))
771        allocate(zdtfic_omp(klon,llm))
772        allocate(zdqfic_omp(klon,llm,nqtot))
773        allocate(zdpsrf_omp(klon))
774        allocate(flxwfi_omp(klon,llm))
775
776        if (flag_moyzon) call moyzon_init_omp(klon,llm,nqtot)
777
778        first_omp=.false.
779      endif
780       
781           
782      klon=klon_omp
783      offset=klon_omp_begin-1
784     
785      do l=1,llm+1
786        do i=1,klon
787          zplev_omp(i,l)=zplev(offset+i,l)
788        enddo
789      enddo
790         
791       do l=1,llm
792        do i=1,klon 
793          zplay_omp(i,l)=zplay(offset+i,l)
794        enddo
795      enddo
796       
797       do l=1,llm
798        do i=1,klon 
799          zpk_omp(i,l)=zpk(offset+i,l)
800        enddo
801      enddo
802       
803      do l=1,llm
804        do i=1,klon
805          zphi_omp(i,l)=zphi(offset+i,l)
806        enddo
807      enddo
808       
809      do i=1,klon
810        zphis_omp(i)=zphis(offset+i)
811      enddo
812     
813       
814      do l=1,llm
815        presnivs_omp(l)=presnivs(l)
816      enddo
817       
818      do l=1,llm
819        do i=1,klon
820          zufi_omp(i,l)=zufi(offset+i,l)
821        enddo
822      enddo
823       
824      do l=1,llm
825        do i=1,klon
826          zvfi_omp(i,l)=zvfi(offset+i,l)
827        enddo
828      enddo
829       
830      do l=1,llm
831        do i=1,klon
832          zrfi_omp(i,l)=zrfi(offset+i,l)
833        enddo
834      enddo
835       
836       
837      do l=1,llm
838        do i=1,klon
839          ztfi_omp(i,l)=ztfi(offset+i,l)
840        enddo
841      enddo
842       
843      do iq=1,nqtot
844        do l=1,llm
845          do i=1,klon
846            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
847          enddo
848        enddo
849      enddo
850       
851      do l=1,llm
852        do i=1,klon
853          zdufi_omp(i,l)=zdufi(offset+i,l)
854        enddo
855      enddo
856       
857      do l=1,llm
858        do i=1,klon
859          zdvfi_omp(i,l)=zdvfi(offset+i,l)
860        enddo
861      enddo
862       
863      do l=1,llm
864        do i=1,klon
865          zdtfi_omp(i,l)=zdtfi(offset+i,l)
866        enddo
867      enddo
868       
869      do iq=1,nqtot
870        do l=1,llm
871          do i=1,klon
872            zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
873          enddo
874        enddo
875      enddo
876       
877      do i=1,klon
878        zdpsrf_omp(i)=zdpsrf(offset+i)
879      enddo
880
881      do l=1,llm
882        do i=1,klon
883          flxwfi_omp(i,l)=flxwfi(offset+i,l)
884        enddo
885      enddo
886
887      if (flag_moyzon) then
888       do l=1,llm+1
889        do i=1,klon
890          zplevbar(i,l)=zplevbar_mpi(offset+i,l)
891        enddo
892       enddo
893         
894       do l=1,llm
895        do i=1,klon 
896          zplaybar(i,l)=zplaybar_mpi(offset+i,l)
897        enddo
898       enddo
899       
900       do l=1,llm
901        do i=1,klon
902          zphibar(i,l)=zphibar_mpi(offset+i,l)
903        enddo
904       enddo
905               
906      do i=1,klon
907        zphisbar(i)=zphisbar_mpi(offset+i)
908      enddo
909     
910      do l=1,llm
911        do i=1,klon
912          ztfibar(i,l)=ztfibar_mpi(offset+i,l)
913        enddo
914      enddo
915       
916      do iq=1,nqtot
917        do l=1,llm
918          do i=1,klon
919            zqfibar(i,l,iq)=zqfibar_mpi(offset+i,l,iq)
920          enddo
921        enddo
922       enddo
923      endif
924
925     
926c$OMP BARRIER
927
928!$OMP MASTER
929!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
930!$OMP END MASTER
931      zdt_split=dtphys/nsplit_phys
932      zdufic_omp(:,:)=0.
933      zdvfic_omp(:,:)=0.
934      zdtfic_omp(:,:)=0.
935      zdqfic_omp(:,:,:)=0.
936
937#ifdef CPP_PHYS
938      do isplit=1,nsplit_phys
939
940         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
941         debut_split=debut.and.isplit==1
942         lafin_split=lafin.and.isplit==nsplit_phys
943
944        CALL call_physiq(klon,llm,nqtot,tname,
945     &                   debut_split,lafin_split,
946     &                   jD_cur,jH_cur_split,zdt_split,
947     &                   zplev_omp,zplay_omp,
948     &                   zpk_omp,zphi_omp,zphis_omp,
949     &                   presnivs_omp,
950     &                   zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
951     &                   flxwfi_omp,pducov,
952     &                   zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
953     &                   zdpsrf_omp)
954
955!      if (planet_type=="earth") then
956!        CALL physiq (klon,
957!     .             llm,
958!     .             debut_split,
959!     .             lafin_split,
960!     .             jD_cur,
961!     .             jH_cur_split,
962!     .             zdt_split,
963!     .             zplev_omp,
964!     .             zplay_omp,
965!     .             zphi_omp,
966!     .             zphis_omp,
967!     .             presnivs_omp,
968!     .             zufi_omp,
969!     .             zvfi_omp,
970!     .             ztfi_omp,
971!     .             zqfi_omp,
972!     .             flxwfi_omp,
973!     .             zdufi_omp,
974!     .             zdvfi_omp,
975!     .             zdtfi_omp,
976!     .             zdqfi_omp,
977!     .             zdpsrf_omp,
978!     .             pducov)
979!
980!      else if ( planet_type=="generic" ) then
981!
982!      CALL physiq (klon,     !! ngrid
983!     .             llm,            !! nlayer
984!     .             nqtot,          !! nq
985!     .             tname,          !! tracer names from dynamical core (given in infotrac)
986!     .             debut_split,    !! firstcall
987!     .             lafin_split,    !! lastcall
988!     .             jD_cur,         !! pday. see leapfrog_p
989!     .             jH_cur_split,   !! ptime "fraction of day"
990!     .             zdt_split,      !! ptimestep
991!     .             zplev_omp,  !! pplev
992!     .             zplay_omp,  !! pplay
993!     .             zphi_omp,   !! pphi
994!     .             zufi_omp,   !! pu
995!     .             zvfi_omp,   !! pv
996!     .             ztfi_omp,   !! pt
997!     .             zqfi_omp,   !! pq
998!     .             flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
999!     .             zdufi_omp,  !! pdu
1000!     .             zdvfi_omp,  !! pdv
1001!     .             zdtfi_omp,  !! pdt
1002!     .             zdqfi_omp,  !! pdq
1003!     .             zdpsrf_omp, !! pdpsrf
1004!     .             tracerdyn)      !! tracerdyn <-- utilite ???
1005!
1006!      else if ( planet_type=="mars" ) then
1007!
1008!        CALL physiq (klon,       ! ngrid
1009!     .             llm,          ! nlayer
1010!     .             nqtot,        ! nq
1011!     .             debut_split,  ! firstcall
1012!     .             lafin_split,  ! lastcall
1013!     .             jD_cur,       ! pday
1014!     .             jH_cur_split, ! ptime
1015!     .             zdt_split,    ! ptimestep
1016!     .             zplev_omp,    ! pplev
1017!     .             zplay_omp,    ! pplay
1018!     .             zphi_omp,     ! pphi
1019!     .             zufi_omp,     ! pu
1020!     .             zvfi_omp,     ! pv
1021!     .             ztfi_omp,     ! pt
1022!     .             zqfi_omp,     ! pq
1023!     .             flxwfi_omp,   ! pw
1024!     .             zdufi_omp,    ! pdu
1025!     .             zdvfi_omp,    ! pdv
1026!     .             zdtfi_omp,    ! pdt
1027!     .             zdqfi_omp,    ! pdq
1028!     .             zdpsrf_omp,   ! pdpsrf
1029!     .             tracerdyn)    ! tracerdyn (somewhat obsolete)
1030!
1031!      else if ((planet_type=="titan").or.(planet_type=="venus")) then
1032!
1033!        CALL physiq (klon,
1034!     .             llm,
1035!     .             nqtot,
1036!     .             debut_split,
1037!     .             lafin_split,
1038!     .             jD_cur,
1039!     .             jH_cur_split,
1040!     .             zdt_split,
1041!     .             zplev_omp,
1042!     .             zplay_omp,
1043!     .             zpk_omp,
1044!     .             zphi_omp,
1045!     .             zphis_omp,
1046!     .             presnivs_omp,
1047!     .             zufi_omp,
1048!     .             zvfi_omp,
1049!     .             ztfi_omp,
1050!     .             zqfi_omp,
1051!     .             flxwfi_omp,
1052!     .             zdufi_omp,
1053!     .             zdvfi_omp,
1054!     .             zdtfi_omp,
1055!     .             zdqfi_omp,
1056!     .             zdpsrf_omp)
1057!
1058!      else ! unknown "planet_type"
1059!
1060!        write(lunout,*) "calfis_p: error, unknown planet_type: ",
1061!     &                  trim(planet_type)
1062!        stop
1063!
1064!      endif ! planet_type
1065         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
1066         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
1067         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
1068         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
1069
1070         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
1071         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
1072         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
1073         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
1074
1075      enddo
1076
1077#endif
1078! of #ifdef CPP_PHYS
1079
1080      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
1081      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
1082      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
1083      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
1084
1085c$OMP BARRIER
1086
1087      do l=1,llm+1
1088        do i=1,klon
1089          zplev(offset+i,l)=zplev_omp(i,l)
1090        enddo
1091      enddo
1092         
1093       do l=1,llm
1094        do i=1,klon 
1095          zplay(offset+i,l)=zplay_omp(i,l)
1096        enddo
1097      enddo
1098       
1099      do l=1,llm
1100        do i=1,klon
1101          zphi(offset+i,l)=zphi_omp(i,l)
1102        enddo
1103      enddo
1104       
1105
1106      do i=1,klon
1107        zphis(offset+i)=zphis_omp(i)
1108      enddo
1109     
1110       
1111      do l=1,llm
1112        presnivs(l)=presnivs_omp(l)
1113      enddo
1114       
1115      do l=1,llm
1116        do i=1,klon
1117          zufi(offset+i,l)=zufi_omp(i,l)
1118        enddo
1119      enddo
1120       
1121      do l=1,llm
1122        do i=1,klon
1123          zvfi(offset+i,l)=zvfi_omp(i,l)
1124        enddo
1125      enddo
1126       
1127      do l=1,llm
1128        do i=1,klon
1129          ztfi(offset+i,l)=ztfi_omp(i,l)
1130        enddo
1131      enddo
1132       
1133      do iq=1,nqtot
1134        do l=1,llm
1135          do i=1,klon
1136            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
1137          enddo
1138        enddo
1139      enddo
1140       
1141      do l=1,llm
1142        do i=1,klon
1143          zdufi(offset+i,l)=zdufi_omp(i,l)
1144        enddo
1145      enddo
1146       
1147      do l=1,llm
1148        do i=1,klon
1149          zdvfi(offset+i,l)=zdvfi_omp(i,l)
1150        enddo
1151      enddo
1152       
1153      do l=1,llm
1154        do i=1,klon
1155          zdtfi(offset+i,l)=zdtfi_omp(i,l)
1156        enddo
1157      enddo
1158       
1159      do iq=1,nqtot
1160        do l=1,llm
1161          do i=1,klon
1162            zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
1163          enddo
1164        enddo
1165      enddo
1166       
1167      do i=1,klon
1168        zdpsrf(offset+i)=zdpsrf_omp(i)
1169      enddo
1170     
1171
1172      klon=klon_mpi
1173500   CONTINUE
1174c$OMP BARRIER
1175
1176c$OMP MASTER
1177      call stop_timer(timer_physic)
1178c$OMP END MASTER
1179
1180      IF (using_mpi) THEN
1181           
1182      if (MPI_rank>0) then
1183
1184c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
1185       DO l=1,llm     
1186        du_send(1:iim,l)=zdufi(1:iim,l)
1187        dv_send(1:iim,l)=zdvfi(1:iim,l)
1188       ENDDO
1189c$OMP END DO NOWAIT       
1190
1191c$OMP BARRIER
1192#ifdef CPP_MPI
1193c$OMP MASTER
1194!$OMP CRITICAL (MPI)
1195        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
1196     &                   COMM_LMDZ,Req(1),ierr)
1197        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
1198     &                  COMM_LMDZ,Req(2),ierr)
1199!$OMP END CRITICAL (MPI)
1200c$OMP END MASTER
1201#endif
1202c$OMP BARRIER
1203     
1204      endif
1205   
1206      if (MPI_rank<MPI_Size-1) then
1207c$OMP BARRIER
1208#ifdef CPP_MPI
1209c$OMP MASTER     
1210!$OMP CRITICAL (MPI)
1211        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
1212     &                 COMM_LMDZ,Req(3),ierr)
1213        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
1214     &                 COMM_LMDZ,Req(4),ierr)
1215!$OMP END CRITICAL (MPI)
1216c$OMP END MASTER
1217#endif
1218      endif
1219
1220c$OMP BARRIER
1221
1222
1223#ifdef CPP_MPI
1224c$OMP MASTER   
1225!$OMP CRITICAL (MPI)
1226      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
1227        call MPI_WAITALL(4,Req(1),Status,ierr)
1228      else if (MPI_rank>0) then
1229        call MPI_WAITALL(2,Req(1),Status,ierr)
1230      else if (MPI_rank <MPI_Size-1) then
1231        call MPI_WAITALL(2,Req(3),Status,ierr)
1232      endif
1233!$OMP END CRITICAL (MPI)
1234c$OMP END MASTER
1235#endif
1236
1237c$OMP BARRIER     
1238
1239      ENDIF ! using_mpi
1240     
1241     
1242c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1243      DO l=1,llm
1244           
1245        zdufi2(1:klon,l)=zdufi(1:klon,l)
1246        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
1247           
1248        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
1249        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
1250 
1251         pdhfi(:,jj_begin,l)=0
1252         pdqfi(:,jj_begin,l,:)=0
1253         pdufi(:,jj_begin,l)=0
1254         pdvfi(:,jj_begin,l)=0
1255         
1256         if (.not. is_south_pole_dyn) then
1257           pdhfi(:,jj_end,l)=0
1258           pdqfi(:,jj_end,l,:)=0
1259           pdufi(:,jj_end,l)=0
1260           pdvfi(:,jj_end,l)=0
1261         endif
1262     
1263       ENDDO
1264c$OMP END DO NOWAIT
1265
1266c$OMP MASTER
1267       pdpsfi(:,jj_begin)=0   
1268       if (.not. is_south_pole_dyn) then
1269         pdpsfi(:,jj_end)=0
1270       endif
1271c$OMP END MASTER
1272c-----------------------------------------------------------------------
1273c   transformation des tendances physiques en tendances dynamiques:
1274c   ---------------------------------------------------------------
1275
1276c  tendance sur la pression :
1277c  -----------------------------------
1278      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
1279c
1280c   62. enthalpie potentielle
1281c   ---------------------
1282
1283      kstart=1
1284      kend=klon
1285
1286      if (is_north_pole_dyn) kstart=2
1287      if (is_south_pole_dyn)  kend=klon-1
1288     
1289! ADAPTATION GCM POUR CP(T)
1290      call t2tpot_p(klon,llm,ztfi,zteta,zpk)
1291
1292
1293c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1294      DO l=1,llm
1295!CDIR ON_ADB(index_i)
1296!CDIR ON_ADB(index_j)
1297!cdir NODEP
1298        do ig0=kstart,kend
1299          i=index_i(ig0)
1300          j=index_j(ig0)
1301!          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1302          pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1303!          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
1304          if (i==1) then
1305            pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1306          endif
1307         enddo         
1308
1309        if (is_north_pole_dyn) then
1310            DO i=1,iip1
1311!              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
1312              pdhfi(i,1,l)    = (zteta(1,l) - pteta(i,1,l))/dtphys
1313            enddo
1314        endif
1315       
1316        if (is_south_pole_dyn) then
1317            DO i=1,iip1
1318!              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
1319              pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
1320            ENDDO
1321        endif
1322      ENDDO
1323c$OMP END DO NOWAIT
1324     
1325c   62. humidite specifique
1326c   ---------------------
1327! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1328!      DO iq=1,nqtot
1329!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1330!         DO l=1,llm
1331!!!cdir NODEP
1332!           do ig0=kstart,kend
1333!             i=index_i(ig0)
1334!             j=index_j(ig0)
1335!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1336!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1337!           enddo
1338!           
1339!           if (is_north_pole_dyn) then
1340!             do i=1,iip1
1341!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
1342!             enddo
1343!           endif
1344!           
1345!           if (is_south_pole_dyn) then
1346!             do i=1,iip1
1347!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1348!             enddo
1349!           endif
1350!         ENDDO
1351!c$OMP END DO NOWAIT
1352!      ENDDO
1353
1354c   63. traceurs (tous en intensifs)
1355c   ------------
1356C     initialisation des tendances
1357
1358c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1359      DO l=1,llm
1360        pdqfi(:,:,l,:)=0.
1361      ENDDO
1362c$OMP END DO NOWAIT     
1363
1364C
1365!cdir NODEP
1366      DO iq=1,nqtot
1367c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1368         DO l=1,llm
1369!CDIR ON_ADB(index_i)
1370!CDIR ON_ADB(index_j)
1371!cdir NODEP           
1372             DO ig0=kstart,kend
1373              i=index_i(ig0)
1374              j=index_j(ig0)
1375              pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1376              if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1377            ENDDO
1378           
1379            IF (is_north_pole_dyn) then
1380              DO i=1,iip1
1381                pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
1382              ENDDO
1383            ENDIF
1384           
1385            IF (is_south_pole_dyn) then
1386              DO i=1,iip1
1387                pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1388              ENDDO
1389            ENDIF
1390           
1391         ENDDO
1392c$OMP END DO NOWAIT     
1393      ENDDO
1394     
1395c   65. champ u:
1396c   ------------
1397c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1398      DO l=1,llm
1399!CDIR ON_ADB(index_i)
1400!CDIR ON_ADB(index_j)
1401!cdir NODEP
1402         do ig0=kstart,kend
1403           i=index_i(ig0)
1404           j=index_j(ig0)
1405           
1406           if (i/=iim) then
1407             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1408           endif
1409           
1410           if (i==1) then
1411              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
1412     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
1413             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1414           endif
1415         
1416         enddo
1417         
1418         if (is_north_pole_dyn) then
1419           DO i=1,iip1
1420            pdufi(i,1,l)    = 0.
1421           ENDDO
1422         endif
1423         
1424         if (is_south_pole_dyn) then
1425           DO i=1,iip1
1426            pdufi(i,jjp1,l) = 0.
1427           ENDDO
1428         endif
1429         
1430      ENDDO
1431c$OMP END DO NOWAIT
1432
1433c   67. champ v:
1434c   ------------
1435
1436      kstart=1
1437      kend=klon
1438
1439      if (is_north_pole_dyn) kstart=2
1440      if (is_south_pole_dyn)  kend=klon-1-iim
1441     
1442c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1443      DO l=1,llm
1444!CDIR ON_ADB(index_i)
1445!CDIR ON_ADB(index_j)
1446!cdir NODEP
1447        do ig0=kstart,kend
1448           i=index_i(ig0)
1449           j=index_j(ig0)
1450           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
1451           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1452     $                                      zdvfi2(ig0+iim,l))
1453     $                                    *cv(i,j)
1454        enddo
1455         
1456      ENDDO
1457c$OMP END DO NOWAIT
1458
1459
1460c   68. champ v pres des poles:
1461c   ---------------------------
1462c      v = U * cos(long) + V * SIN(long)
1463
1464      if (is_north_pole_dyn) then
1465
1466c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1467        DO l=1,llm
1468
1469          DO i=1,iim
1470            pdvfi(i,1,l)=
1471     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
1472       
1473            pdvfi(i,1,l)=
1474     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
1475          ENDDO
1476
1477          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
1478
1479        ENDDO
1480c$OMP END DO NOWAIT
1481
1482      endif   
1483     
1484      if (is_south_pole_dyn) then
1485
1486c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1487         DO l=1,llm
1488 
1489           DO i=1,iim
1490              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1491     $        +zdvfi(klon,l)*SIN(rlonv(i))
1492
1493              pdvfi(i,jjm,l)=
1494     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
1495           ENDDO
1496
1497           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1498
1499        ENDDO
1500c$OMP END DO NOWAIT
1501     
1502      endif
1503c-----------------------------------------------------------------------
1504
1505700   CONTINUE
1506 
1507      firstcal = .FALSE.
1508
1509#else
1510      write(lunout,*)
1511     & "calfis_p: for now can only work with parallel physics"
1512      stop
1513#endif
1514! of #ifdef CPP_PHYS
1515#endif
1516! of #ifdef CPP_PARA
1517      END
Note: See TracBrowser for help on using the repository browser.