source: trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/calfis_p.F @ 1459

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

Common dynamics: A couple of bug fixes

  • calfis[_p].F array boundaries must be explicitely specified when underlying arrays are of different sizes.
  • advect_new_p.F : missing initializations of intermediate variables at topmost layer.

EM

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