source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F @ 245

Last change on this file since 245 was 245, checked in by ymipsl, 10 years ago
  • One call for initialize physics from dynamico
  • mpi_root renamed into mpi_master due to conflict with an existaing symbol from the mpi library

==> mpi_root => mpi_master, is_mpi_root => is_mpi_master, is_omp_root => is_omp_master

YM
File size: 39.9 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
32      USE mod_interface_dyn_phys
33      USE mod_const_mpi, ONLY : COMM_LMDZ
34!      USE IOPHY
35#endif
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 infotrac, ONLY: nqtot, niadv, tname
41      USE control_mod, ONLY: planet_type, nsplit_phys
42      USE cpdet_mod, only: tpot2t_p, t2tpot_p
43
44! used only for zonal averages
45      USE moyzon_mod
46
47      IMPLICIT NONE
48c=======================================================================
49c
50c   1. rearrangement des tableaux et transformation
51c      variables dynamiques  >  variables physiques
52c   2. calcul des termes physiques
53c   3. retransformation des tendances physiques en tendances dynamiques
54c
55c   remarques:
56c   ----------
57c
58c    - les vents sont donnes dans la physique par leurs composantes 
59c      naturelles.
60c    - la variable thermodynamique de la physique est une variable
61c      intensive :   T 
62c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
63c    - les deux seules variables dependant de la geometrie necessaires
64c      pour la physique sont la latitude pour le rayonnement et 
65c      l'aire de la maille quand on veut integrer une grandeur
66c      horizontalement.
67c    - les points de la physique sont les points scalaires de la
68c      la dynamique; numerotation:
69c          1 pour le pole nord
70c          (jjm-1)*iim pour l'interieur du domaine
71c          ngridmx pour le pole sud
72c      ---> ngridmx=2+(jjm-1)*iim
73c
74c     Input :
75c     -------
76c       ecritphy        frequence d'ecriture (en jours)de histphy
77c       pucov           covariant zonal velocity
78c       pvcov           covariant meridional velocity
79c       pteta           potential temperature
80c       pps             surface pressure
81c       pmasse          masse d'air dans chaque maille
82c       pts             surface temperature  (K)
83c       callrad         clef d'appel au rayonnement
84c
85c    Output :
86c    --------
87c        pdufi          tendency for the natural zonal velocity (ms-1)
88c        pdvfi          tendency for the natural meridional velocity
89c        pdhfi          tendency for the potential temperature (K/s)
90c        pdtsfi         tendency for the surface temperature
91c
92c        pdtrad         radiative tendencies  \  both input
93c        pfluxrad       radiative fluxes      /  and output
94c
95c=======================================================================
96c
97c-----------------------------------------------------------------------
98c
99c    0.  Declarations :
100c    ------------------
101
102#include "dimensions.h"
103#include "paramet.h"
104#include "temps.h"
105#include "logic.h"
106
107      INTEGER ngridmx
108      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
109
110#include "comconst.h"
111#include "comvert.h"
112#include "comgeom2.h"
113#include "iniprint.h"
114#ifdef CPP_MPI
115      include 'mpif.h'
116#endif
117c    Arguments :
118c    -----------
119      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
120      REAL,INTENT(IN) :: jD_cur, jH_cur
121      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
122      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
123      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
124      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
125      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
126      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
127      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
128
129      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
130      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
131      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
132! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
133! qui lui meme ne sert a rien dans la routine telle qu'elle est
134! ecrite, et que j'ai donc commente....
135      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
136      ! NB: pdq is only used to compute pcvgq which is in fact not used...
137
138      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
139      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
140      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
141      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm)  ! Vertical mass flux on dynamics grid
142
143      ! tendencies (in */s) from the physics
144      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
145      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
146      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
147      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
148      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
149
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(:) = g*rad*rad/(g*rad-phimoy(:))-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,
898c#ifdef INCA
899     .             flxwfi_omp,
900c#endif
901     .             zdufi_omp,
902     .             zdvfi_omp,
903     .             zdtfi_omp,
904     .             zdqfi_omp,
905     .             zdpsrf_omp,
906     .             pducov)
907
908      else if ( planet_type=="generic" ) then
909
910      CALL physiq (klon,     !! ngrid
911     .             llm,            !! nlayer
912     .             nqtot,          !! nq
913     .             tname,          !! tracer names from dynamical core (given in infotrac)
914     .             debut_split,    !! firstcall
915     .             lafin_split,    !! lastcall
916     .             jD_cur,         !! pday. see leapfrog_p
917     .             jH_cur_split,   !! ptime "fraction of day"
918     .             zdt_split,      !! ptimestep
919     .             zplev_omp,  !! pplev
920     .             zplay_omp,  !! pplay
921     .             zphi_omp,   !! pphi
922     .             zufi_omp,   !! pu
923     .             zvfi_omp,   !! pv
924     .             ztfi_omp,   !! pt
925     .             zqfi_omp,   !! pq
926     .             flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
927     .             zdufi_omp,  !! pdu
928     .             zdvfi_omp,  !! pdv
929     .             zdtfi_omp,  !! pdt
930     .             zdqfi_omp,  !! pdq
931     .             zdpsrf_omp, !! pdpsrf
932     .             tracerdyn)      !! tracerdyn <-- utilite ???
933
934      else if ( planet_type=="mars" ) then
935
936        CALL physiq (klon,       ! ngrid
937     .             llm,          ! nlayer
938     .             nqtot,        ! nq
939     .             debut_split,  ! firstcall
940     .             lafin_split,  ! lastcall
941     .             jD_cur,       ! pday
942     .             jH_cur_split, ! ptime
943     .             zdt_split,    ! ptimestep
944     .             zplev_omp,    ! pplev
945     .             zplay_omp,    ! pplay
946     .             zphi_omp,     ! pphi
947     .             zufi_omp,     ! pu
948     .             zvfi_omp,     ! pv
949     .             ztfi_omp,     ! pt
950     .             zqfi_omp,     ! pq
951     .             flxwfi_omp,   ! pw
952     .             zdufi_omp,    ! pdu
953     .             zdvfi_omp,    ! pdv
954     .             zdtfi_omp,    ! pdt
955     .             zdqfi_omp,    ! pdq
956     .             zdpsrf_omp,   ! pdpsrf
957     .             tracerdyn)    ! tracerdyn (somewhat obsolete)
958
959      else if ((planet_type=="titan").or.(planet_type=="venus")) then
960
961        CALL physiq (klon,
962     .             llm,
963     .             nqtot,
964     .             debut_split,
965     .             lafin_split,
966     .             jD_cur,
967     .             jH_cur_split,
968     .             zdt_split,
969     .             zplev_omp,
970     .             zplay_omp,
971     .             zpk_omp,
972     .             zphi_omp,
973     .             zphis_omp,
974     .             presnivs_omp,
975     .             zufi_omp,
976     .             zvfi_omp,
977     .             ztfi_omp,
978     .             zqfi_omp,
979     .             flxwfi_omp,
980     .             zdufi_omp,
981     .             zdvfi_omp,
982     .             zdtfi_omp,
983     .             zdqfi_omp,
984     .             zdpsrf_omp)
985
986      else ! unknown "planet_type"
987
988        write(lunout,*) "calfis_p: error, unknown planet_type: ",
989     &                  trim(planet_type)
990        stop
991
992      endif ! planet_type
993         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
994         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
995         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
996         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
997
998         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
999         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
1000         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
1001         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
1002
1003      enddo
1004
1005#endif
1006! of #ifdef CPP_PHYS
1007
1008      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
1009      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
1010      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
1011      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
1012
1013c$OMP BARRIER
1014
1015      do l=1,llm+1
1016        do i=1,klon
1017          zplev(offset+i,l)=zplev_omp(i,l)
1018        enddo 
1019      enddo
1020         
1021       do l=1,llm
1022        do i=1,klon 
1023          zplay(offset+i,l)=zplay_omp(i,l)
1024        enddo 
1025      enddo
1026       
1027      do l=1,llm
1028        do i=1,klon
1029          zphi(offset+i,l)=zphi_omp(i,l)
1030        enddo 
1031      enddo
1032       
1033
1034      do i=1,klon
1035        zphis(offset+i)=zphis_omp(i)
1036      enddo 
1037     
1038       
1039      do l=1,llm
1040        presnivs(l)=presnivs_omp(l)
1041      enddo 
1042       
1043      do l=1,llm
1044        do i=1,klon
1045          zufi(offset+i,l)=zufi_omp(i,l)
1046        enddo 
1047      enddo
1048       
1049      do l=1,llm
1050        do i=1,klon
1051          zvfi(offset+i,l)=zvfi_omp(i,l)
1052        enddo 
1053      enddo
1054       
1055      do l=1,llm
1056        do i=1,klon
1057          ztfi(offset+i,l)=ztfi_omp(i,l)
1058        enddo 
1059      enddo
1060       
1061      do iq=1,nqtot
1062        do l=1,llm
1063          do i=1,klon
1064            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
1065          enddo
1066        enddo 
1067      enddo
1068       
1069      do l=1,llm
1070        do i=1,klon
1071          zdufi(offset+i,l)=zdufi_omp(i,l)
1072        enddo 
1073      enddo
1074       
1075      do l=1,llm
1076        do i=1,klon
1077          zdvfi(offset+i,l)=zdvfi_omp(i,l)
1078        enddo 
1079      enddo
1080       
1081      do l=1,llm
1082        do i=1,klon
1083          zdtfi(offset+i,l)=zdtfi_omp(i,l)
1084        enddo 
1085      enddo
1086       
1087      do iq=1,nqtot
1088        do l=1,llm
1089          do i=1,klon
1090            zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
1091          enddo 
1092        enddo
1093      enddo
1094       
1095      do i=1,klon
1096        zdpsrf(offset+i)=zdpsrf_omp(i)
1097      enddo 
1098     
1099
1100      klon=klon_mpi
1101500   CONTINUE
1102c$OMP BARRIER
1103
1104c$OMP MASTER
1105      call stop_timer(timer_physic)
1106c$OMP END MASTER
1107
1108      IF (using_mpi) THEN
1109           
1110      if (MPI_rank>0) then
1111
1112c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
1113       DO l=1,llm     
1114        du_send(1:iim,l)=zdufi(1:iim,l)
1115        dv_send(1:iim,l)=zdvfi(1:iim,l)
1116       ENDDO
1117c$OMP END DO NOWAIT       
1118
1119c$OMP BARRIER
1120#ifdef CPP_MPI 
1121c$OMP MASTER
1122!$OMP CRITICAL (MPI)
1123        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
1124     &                   COMM_LMDZ,Req(1),ierr)
1125        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
1126     &                  COMM_LMDZ,Req(2),ierr)
1127!$OMP END CRITICAL (MPI)
1128c$OMP END MASTER
1129#endif
1130c$OMP BARRIER
1131     
1132      endif
1133   
1134      if (MPI_rank<MPI_Size-1) then
1135c$OMP BARRIER
1136#ifdef CPP_MPI 
1137c$OMP MASTER     
1138!$OMP CRITICAL (MPI)
1139        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
1140     &                 COMM_LMDZ,Req(3),ierr)
1141        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
1142     &                 COMM_LMDZ,Req(4),ierr)
1143!$OMP END CRITICAL (MPI)
1144c$OMP END MASTER
1145#endif
1146      endif
1147
1148c$OMP BARRIER
1149
1150
1151#ifdef CPP_MPI 
1152c$OMP MASTER   
1153!$OMP CRITICAL (MPI)
1154      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
1155        call MPI_WAITALL(4,Req(1),Status,ierr)
1156      else if (MPI_rank>0) then
1157        call MPI_WAITALL(2,Req(1),Status,ierr)
1158      else if (MPI_rank <MPI_Size-1) then
1159        call MPI_WAITALL(2,Req(3),Status,ierr)
1160      endif
1161!$OMP END CRITICAL (MPI)
1162c$OMP END MASTER
1163#endif
1164
1165c$OMP BARRIER     
1166
1167      ENDIF ! using_mpi
1168     
1169     
1170c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1171      DO l=1,llm
1172           
1173        zdufi2(1:klon,l)=zdufi(1:klon,l)
1174        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
1175           
1176        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
1177        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l) 
1178 
1179         pdhfi(:,jj_begin,l)=0
1180         pdqfi(:,jj_begin,l,:)=0
1181         pdufi(:,jj_begin,l)=0
1182         pdvfi(:,jj_begin,l)=0
1183         
1184         if (.not. is_south_pole) then
1185           pdhfi(:,jj_end,l)=0
1186           pdqfi(:,jj_end,l,:)=0
1187           pdufi(:,jj_end,l)=0
1188           pdvfi(:,jj_end,l)=0
1189         endif
1190     
1191       ENDDO 
1192c$OMP END DO NOWAIT
1193
1194c$OMP MASTER
1195       pdpsfi(:,jj_begin)=0   
1196       if (.not. is_south_pole) then
1197         pdpsfi(:,jj_end)=0
1198       endif
1199c$OMP END MASTER
1200c-----------------------------------------------------------------------
1201c   transformation des tendances physiques en tendances dynamiques:
1202c   ---------------------------------------------------------------
1203
1204c  tendance sur la pression :
1205c  -----------------------------------
1206      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
1207c
1208c   62. enthalpie potentielle
1209c   ---------------------
1210
1211      kstart=1
1212      kend=klon
1213
1214      if (is_north_pole) kstart=2
1215      if (is_south_pole)  kend=klon-1
1216     
1217! ADAPTATION GCM POUR CP(T)
1218      call t2tpot_p(klon,llm,ztfi,zteta,zpk)
1219
1220
1221c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1222      DO l=1,llm
1223!CDIR ON_ADB(index_i)
1224!CDIR ON_ADB(index_j)
1225!cdir NODEP
1226        do ig0=kstart,kend
1227          i=index_i(ig0)
1228          j=index_j(ig0)
1229!          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1230          pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1231!          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
1232          if (i==1) then
1233            pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1234          endif
1235         enddo         
1236
1237        if (is_north_pole) then
1238            DO i=1,iip1
1239!              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
1240              pdhfi(i,1,l)    = (zteta(1,l) - pteta(i,1,l))/dtphys
1241            enddo
1242        endif
1243       
1244        if (is_south_pole) then
1245            DO i=1,iip1
1246!              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
1247              pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
1248            ENDDO
1249        endif
1250      ENDDO
1251c$OMP END DO NOWAIT
1252     
1253c   62. humidite specifique
1254c   ---------------------
1255! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1256!      DO iq=1,nqtot
1257!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1258!         DO l=1,llm
1259!!!cdir NODEP
1260!           do ig0=kstart,kend
1261!             i=index_i(ig0)
1262!             j=index_j(ig0)
1263!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1264!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1265!           enddo
1266!           
1267!           if (is_north_pole) then
1268!             do i=1,iip1
1269!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
1270!             enddo
1271!           endif
1272!           
1273!           if (is_south_pole) then
1274!             do i=1,iip1
1275!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1276!             enddo
1277!           endif
1278!         ENDDO
1279!c$OMP END DO NOWAIT
1280!      ENDDO
1281
1282c   63. traceurs (tous en intensifs)
1283c   ------------
1284C     initialisation des tendances
1285
1286c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1287      DO l=1,llm
1288        pdqfi(:,:,l,:)=0.
1289      ENDDO
1290c$OMP END DO NOWAIT     
1291
1292C
1293!cdir NODEP
1294      DO iq=1,nqtot
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              pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1304              if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1305            ENDDO
1306           
1307            IF (is_north_pole) then
1308              DO i=1,iip1
1309                pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
1310              ENDDO
1311            ENDIF
1312           
1313            IF (is_south_pole) then
1314              DO i=1,iip1
1315                pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1316              ENDDO
1317            ENDIF
1318           
1319         ENDDO
1320c$OMP END DO NOWAIT     
1321      ENDDO
1322     
1323c   65. champ u:
1324c   ------------
1325c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1326      DO l=1,llm
1327!CDIR ON_ADB(index_i)
1328!CDIR ON_ADB(index_j)
1329!cdir NODEP
1330         do ig0=kstart,kend
1331           i=index_i(ig0)
1332           j=index_j(ig0)
1333           
1334           if (i/=iim) then
1335             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1336           endif
1337           
1338           if (i==1) then
1339              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
1340     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
1341             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1342           endif
1343         
1344         enddo
1345         
1346         if (is_north_pole) then
1347           DO i=1,iip1
1348            pdufi(i,1,l)    = 0.
1349           ENDDO
1350         endif
1351         
1352         if (is_south_pole) then
1353           DO i=1,iip1
1354            pdufi(i,jjp1,l) = 0.
1355           ENDDO
1356         endif
1357         
1358      ENDDO
1359c$OMP END DO NOWAIT
1360
1361c   67. champ v:
1362c   ------------
1363
1364      kstart=1
1365      kend=klon
1366
1367      if (is_north_pole) kstart=2
1368      if (is_south_pole)  kend=klon-1-iim
1369     
1370c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1371      DO l=1,llm
1372!CDIR ON_ADB(index_i)
1373!CDIR ON_ADB(index_j)
1374!cdir NODEP
1375        do ig0=kstart,kend
1376           i=index_i(ig0)
1377           j=index_j(ig0)
1378           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
1379           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1380     $                                      zdvfi2(ig0+iim,l))
1381     $                                    *cv(i,j)
1382        enddo
1383         
1384      ENDDO
1385c$OMP END DO NOWAIT
1386
1387
1388c   68. champ v pres des poles:
1389c   ---------------------------
1390c      v = U * cos(long) + V * SIN(long)
1391
1392      if (is_north_pole) then
1393
1394c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1395        DO l=1,llm
1396
1397          DO i=1,iim
1398            pdvfi(i,1,l)=
1399     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
1400       
1401            pdvfi(i,1,l)=
1402     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
1403          ENDDO
1404
1405          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
1406
1407        ENDDO
1408c$OMP END DO NOWAIT
1409
1410      endif   
1411     
1412      if (is_south_pole) then
1413
1414c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1415         DO l=1,llm
1416 
1417           DO i=1,iim
1418              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1419     $        +zdvfi(klon,l)*SIN(rlonv(i))
1420
1421              pdvfi(i,jjm,l)=
1422     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
1423           ENDDO
1424
1425           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1426
1427        ENDDO
1428c$OMP END DO NOWAIT
1429     
1430      endif
1431c-----------------------------------------------------------------------
1432
1433700   CONTINUE
1434 
1435      firstcal = .FALSE.
1436
1437#else
1438      write(lunout,*)
1439     & "calfis_p: for now can only work with parallel physics"
1440      stop
1441#endif
1442! of #ifdef CPP_PHYS
1443      RETURN
1444      END
Note: See TracBrowser for help on using the repository browser.