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

Last change on this file since 3897 was 3897, checked in by emillour, 3 months ago

Dynamics-physics interface for Venus PCM:
Remove computation of mean temperature and presure profiles for
Venus physics as:
1 - They are not used
2 - They cause major issues when turning on OpenMP (which will need

further investigation if ever we want to use OpenMP for Titan).

EM

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