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

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

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