source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/calfis_p.F @ 1200

Last change on this file since 1200 was 1142, checked in by yann meurdesoif, 16 years ago

Portage vers Vargas

Y.M + E.M

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.1 KB
RevLine 
[630]1!
2! $Header$
3!
4C
5C
[1114]6      SUBROUTINE calfis_p(lafin,
[630]7     $                  rdayvrai,
8     $                  heure,
9     $                  pucov,
10     $                  pvcov,
11     $                  pteta,
12     $                  pq,
13     $                  pmasse,
14     $                  pps,
15     $                  pp,
16     $                  ppk,
17     $                  pphis,
18     $                  pphi,
19     $                  pducov,
20     $                  pdvcov,
21     $                  pdteta,
22     $                  pdq,
23     $                  flxw,
24     $                  clesphy0,
25     $                  pdufi,
26     $                  pdvfi,
27     $                  pdhfi,
28     $                  pdqfi,
29     $                  pdpsfi)
30c
31c    Auteur :  P. Le Van, F. Hourdin
32c   .........
33      USE dimphy
[774]34      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
[1000]35      USE parallel, ONLY : omp_chunk, using_mpi
[774]36      USE mod_interface_dyn_phys
[630]37      USE Write_Field
38      Use Write_field_p
39      USE Times
[764]40      USE IOPHY
[1114]41      USE infotrac
42
[630]43      IMPLICIT NONE
44c=======================================================================
45c
46c   1. rearrangement des tableaux et transformation
47c      variables dynamiques  >  variables physiques
48c   2. calcul des termes physiques
49c   3. retransformation des tendances physiques en tendances dynamiques
50c
51c   remarques:
52c   ----------
53c
54c    - les vents sont donnes dans la physique par leurs composantes
55c      naturelles.
56c    - la variable thermodynamique de la physique est une variable
57c      intensive :   T
58c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
59c    - les deux seules variables dependant de la geometrie necessaires
60c      pour la physique sont la latitude pour le rayonnement et
61c      l'aire de la maille quand on veut integrer une grandeur
62c      horizontalement.
63c    - les points de la physique sont les points scalaires de la
64c      la dynamique; numerotation:
65c          1 pour le pole nord
66c          (jjm-1)*iim pour l'interieur du domaine
67c          ngridmx pour le pole sud
68c      ---> ngridmx=2+(jjm-1)*iim
69c
70c     Input :
71c     -------
72c       ecritphy        frequence d'ecriture (en jours)de histphy
73c       pucov           covariant zonal velocity
74c       pvcov           covariant meridional velocity
75c       pteta           potential temperature
76c       pps             surface pressure
77c       pmasse          masse d'air dans chaque maille
78c       pts             surface temperature  (K)
79c       callrad         clef d'appel au rayonnement
80c
81c    Output :
82c    --------
83c        pdufi          tendency for the natural zonal velocity (ms-1)
84c        pdvfi          tendency for the natural meridional velocity
85c        pdhfi          tendency for the potential temperature
86c        pdtsfi         tendency for the surface temperature
87c
88c        pdtrad         radiative tendencies  \  both input
89c        pfluxrad       radiative fluxes      /  and output
90c
91c=======================================================================
92c
93c-----------------------------------------------------------------------
94c
95c    0.  Declarations :
96c    ------------------
97
98#include "dimensions.h"
99#include "paramet.h"
100#include "temps.h"
101
[1114]102      INTEGER ngridmx
[630]103      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
104
105#include "comconst.h"
106#include "comvert.h"
107#include "comgeom2.h"
108#include "control.h"
[1000]109#ifdef CPP_MPI
[630]110      include 'mpif.h'
[1000]111#endif
[630]112c    Arguments :
113c    -----------
114      LOGICAL  lafin
115      REAL heure
116
117      REAL pvcov(iip1,jjm,llm)
118      REAL pucov(iip1,jjp1,llm)
119      REAL pteta(iip1,jjp1,llm)
120      REAL pmasse(iip1,jjp1,llm)
[1114]121      REAL pq(iip1,jjp1,llm,nqtot)
[630]122      REAL pphis(iip1,jjp1)
123      REAL pphi(iip1,jjp1,llm)
124c
125      REAL pdvcov(iip1,jjm,llm)
126      REAL pducov(iip1,jjp1,llm)
127      REAL pdteta(iip1,jjp1,llm)
[1114]128      REAL pdq(iip1,jjp1,llm,nqtot)
[630]129c
130      REAL pps(iip1,jjp1)
131      REAL pp(iip1,jjp1,llmp1)
132      REAL ppk(iip1,jjp1,llm)
133c
134      REAL pdvfi(iip1,jjm,llm)
135      REAL pdufi(iip1,jjp1,llm)
136      REAL pdhfi(iip1,jjp1,llm)
[1114]137      REAL pdqfi(iip1,jjp1,llm,nqtot)
[630]138      REAL pdpsfi(iip1,jjp1)
139
140      INTEGER        longcles
141      PARAMETER    ( longcles = 20 )
142      REAL clesphy0( longcles )
143
144
145c    Local variables :
146c    -----------------
147
148      INTEGER i,j,l,ig0,ig,iq,iiq
[764]149      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
150      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
151      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
[630]152c
[764]153      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
154      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
[630]155c
[764]156      REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
157      REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
[630]158c
[960]159c      REAL,ALLOCATABLE,SAVE :: pvervel(:,:)
[630]160c
[764]161      REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
162      REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
163      REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
[985]164      REAL,SAVE,ALLOCATABLE ::  flxwfi(:,:)     ! Flux de masse verticale sur la grille physiq
165
[630]166c
[764]167      REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
168      REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
169      REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
170      REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
171      REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
172      REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
173      REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
174      REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
175      REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
[960]176c      REAL,ALLOCATABLE,SAVE :: pvervel_omp(:,:)
[764]177      REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
178      REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
179      REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
180      REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
181      REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
[985]182      REAL,SAVE,ALLOCATABLE ::  flxwfi_omp(:,:)     ! Flux de masse verticale sur la grille physiq
[764]183
184c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
185c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
[985]186c$OMP+                 zqfi_omp,zdufi_omp,zdvfi_omp,
187c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp)       
[764]188
189      LOGICAL,SAVE :: first_omp=.true.
190c$OMP THREADPRIVATE(first_omp)
191     
[630]192      REAL zsin(iim),zcos(iim),z1(iim)
193      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
194      REAL unskap, pksurcp
[764]195c
196cIM diagnostique PVteta, Amip2
197      INTEGER ntetaSTD
198      PARAMETER(ntetaSTD=3)
199      REAL rtetaSTD(ntetaSTD)
200      DATA rtetaSTD/350., 380., 405./
201      REAL PVteta(klon,ntetaSTD)
202     
[960]203      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
[630]204     
205      REAL SSUM
206
207      LOGICAL firstcal, debut
208      DATA firstcal/.true./
209      SAVE firstcal,debut
[764]210c$OMP THREADPRIVATE(firstcal,debut)
[630]211      REAL rdayvrai
212     
[764]213      REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
[630]214      INTEGER :: ierr
[1000]215#ifdef CPP_MPI
[630]216      INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
[1000]217#else
218      INTEGER,dimension(1,4) :: Status
219#endif
[630]220      INTEGER, dimension(4) :: Req
[764]221      REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
222      integer :: k,kstart,kend
223      INTEGER :: offset 
[630]224c
225c-----------------------------------------------------------------------
226c
227c    1. Initialisations :
228c    --------------------
229c
230
[764]231      klon=klon_mpi
232     
233      PVteta(:,:)=0.
234           
[630]235      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
236         PRINT*,'STOP dans calfis'
237         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
238         PRINT*,'  ngridmx  jjm   iim   '
239         PRINT*,ngridmx,jjm,iim
240         STOP
241      ENDIF
242
243c-----------------------------------------------------------------------
244c   latitude, longitude et aires des mailles pour la physique:
245c   ----------------------------------------------------------
246
247c
248      IF ( firstcal )  THEN
249          debut = .TRUE.
[764]250c$OMP MASTER
251      ALLOCATE(zpsrf(klon))
252      ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
253      ALLOCATE(zphi(klon,llm),zphis(klon))
254      ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
[1114]255      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
[764]256      ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
257      ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
[960]258c      ALLOCATE(pvervel(klon,llm))
[764]259      ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
[1114]260      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
[764]261      ALLOCATE(zdpsrf(klon))
262      ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
[985]263      ALLOCATE(flxwfi(klon,llm))
[764]264c$OMP END MASTER
265c$OMP BARRIER     
[630]266      ELSE
267          debut = .FALSE.
268      ENDIF
269
270c
271c
272c-----------------------------------------------------------------------
273c   40. transformation des variables dynamiques en variables physiques:
274c   ---------------------------------------------------------------
275
276c   41. pressions au sol (en Pascals)
277c   ----------------------------------
278
[764]279c$OMP MASTER
[630]280      call start_timer(timer_physic)
[764]281c$OMP END MASTER
282
283c$OMP MASTER             
[630]284      do ig0=1,klon
[774]285        i=index_i(ig0)
286        j=index_j(ig0)
[630]287        zpsrf(ig0)=pps(i,j)
288      enddo
[764]289c$OMP END MASTER
[630]290
291
292c   42. pression intercouches :
293c
294c   -----------------------------------------------------------------
295c     .... zplev  definis aux (llm +1) interfaces des couches  ....
296c     .... zplay  definis aux (  llm )    milieux des couches  ....
297c   -----------------------------------------------------------------
298
299c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
300c
301       unskap   = 1./ kappa
302c
[961]303c      print *,omp_rank,'klon--->',klon
[764]304c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]305      DO l = 1, llmp1
306        do ig0=1,klon
[774]307          i=index_i(ig0)
308          j=index_j(ig0)
[630]309          zplev( ig0,l ) = pp(i,j,l)
310        enddo
311      ENDDO
[764]312c$OMP END DO NOWAIT
[630]313c
314c
315
316c   43. temperature naturelle (en K) et pressions milieux couches .
317c   ---------------------------------------------------------------
[764]318c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]319      DO l=1,llm
320
321        do ig0=1,klon
[774]322          i=index_i(ig0)
323          j=index_j(ig0)
[630]324          pksurcp        = ppk(i,j,l) / cpp
325          zplay(ig0,l)   = preff * pksurcp ** unskap
326          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
327c          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
328        enddo
329
330      ENDDO
[764]331c$OMP END DO NOWAIT
[630]332
333c   43.bis traceurs
334c   ---------------
335c
336
[1114]337      DO iq=1,nqtot
[630]338         iiq=niadv(iq)
[764]339c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]340         DO l=1,llm
341           do ig0=1,klon
[774]342             i=index_i(ig0)
343             j=index_j(ig0)
[630]344             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
345           enddo
346         ENDDO
[764]347c$OMP END DO NOWAIT     
[630]348      ENDDO
349
350c   convergence dynamique pour les traceurs "EAU"
351
352      DO iq=1,2
[764]353c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]354         DO l=1,llm
355           do ig0=1,klon
[774]356             i=index_i(ig0)
357             j=index_j(ig0)
[630]358c             pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
359           enddo
360         ENDDO
[764]361c$OMP END DO NOWAIT     
[630]362      ENDDO
363
364
365
366c   Geopotentiel calcule par rapport a la surface locale:
367c   -----------------------------------------------------
368
369      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
[764]370
[630]371      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
[1142]372
[764]373c$OMP BARRIER
[630]374
[764]375c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]376      DO l=1,llm
377         DO ig=1,klon
378           zphi(ig,l)=zphi(ig,l)-zphis(ig)
379         ENDDO
380      ENDDO
[960]381c$OMP END DO NOWAIT
382     
[630]383c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
[960]384c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
385c      de masse est calclue dans advtrac_p.F 
[630]386c
[960]387cc$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
388c      DO l=1,llm
389c        do ig0=1,klon
390c           i=index_i(ig0)
391c           j=index_j(ig0)
392c           pvervel(ig0,l) = pw(i,j,l)*g* unsaire(i,j)
393c        enddo
394c       if (is_north_pole) pvervel(1,l)=pw(1,1,l)*g /apoln
395c       if (is_south_pole) pvervel(klon,l)=pw(1,jjp1,l)*g/apols
396c      ENDDO
397cc$OMP END DO NOWAIT
[630]398
399c
400c   45. champ u:
401c   ------------
402
403      kstart=1
404      kend=klon
405     
[774]406      if (is_north_pole) kstart=2
407      if (is_south_pole) kend=klon-1
[630]408     
[764]409c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]410      DO l=1,llm
411        do ig0=kstart,kend
[774]412          i=index_i(ig0)
413          j=index_j(ig0)
[630]414          if (i==1) then
415            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
416     $                         + pucov(1,j,l)/cu(1,j) )
417c            pcvgu(ig0,l)= 0.5*(  pducov(iim,j,l)/cu(iim,j)
418c     $                         + pducov(1,j,l)/cu(1,j) )
419          else
420            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j)
421     $                       + pucov(i,j,l)/cu(i,j) )
422c            pcvgu(ig0,l)= 0.5*(  pducov(i-1,j,l)/cu(i-1,j)
423c     $                        + pducov(i,j,l)/cu(i,j) )
424          endif
425        enddo
426      ENDDO
[764]427c$OMP END DO NOWAIT
[630]428c   46.champ v:
429c   -----------
[764]430c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]431      DO l=1,llm
432        DO ig0=kstart,kend
[774]433          i=index_i(ig0)
434          j=index_j(ig0)
[630]435          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1)
436     $                       + pvcov(i,j,l)/cv(i,j) )
437   
438c          pcvgv(ig0+i,l)= 0.5 * (  pdvcov(i,j-1,l)/cv(i,j-1)
439c     $                           + pdvcov(i,j,l)/cv(i,j) )
440         ENDDO
441      ENDDO
[764]442c$OMP END DO NOWAIT
[630]443
444c   47. champs de vents aux pole nord   
445c   ------------------------------
446c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
447c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
448
[774]449      if (is_north_pole) then
[764]450c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]451        DO l=1,llm
452
453           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
454c           z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
455           DO i=2,iim
456              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
457c              z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
458           ENDDO
459 
460           DO i=1,iim
461              zcos(i)   = COS(rlonv(i))*z1(i)
462c              zcosbis(i)= COS(rlonv(i))*z1bis(i)
463              zsin(i)   = SIN(rlonv(i))*z1(i)
464c              zsinbis(i)= SIN(rlonv(i))*z1bis(i)
465           ENDDO
466 
467           zufi(1,l)  = SSUM(iim,zcos,1)/pi
468c           pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
469           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
470c           pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
471 
472        ENDDO
[764]473c$OMP END DO NOWAIT     
[630]474      endif
475
476
477c   48. champs de vents aux pole sud:
478c   ---------------------------------
479c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
480c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
481
[774]482      if (is_south_pole) then
[764]483c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]484        DO l=1,llm
485 
486         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
487c         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
488           DO i=2,iim
489           z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
490c           z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
491           ENDDO
492 
493           DO i=1,iim
494              zcos(i)    = COS(rlonv(i))*z1(i)
495c              zcosbis(i) = COS(rlonv(i))*z1bis(i)
496              zsin(i)    = SIN(rlonv(i))*z1(i)
497c              zsinbis(i) = SIN(rlonv(i))*z1bis(i)
498           ENDDO
499 
500           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
501c           pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
502           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
503c           pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
504
505        ENDDO
[764]506c$OMP END DO NOWAIT       
[630]507      endif
508
509
[774]510      IF (is_sequential) THEN
[764]511c
512cIM calcul PV a teta=350, 380, 405K
513        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
514     $           ztfi,zplay,zplev,
515     $           ntetaSTD,rtetaSTD,PVteta)
516c
517      ENDIF
[960]518
519c On change de grille, dynamique vers physiq, pour le flux de masse verticale
[630]520      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
521
522c-----------------------------------------------------------------------
523c   Appel de la physique:
524c   ---------------------
525
[764]526cc$OMP  PARALLEL DEFAULT(NONE)
527cc$OMP+ PRIVATE(i,l,offset,iq)
[1114]528cc$OMP+ SHARED(klon_omp_nb,nqtot,klon_omp_begin,
[764]529cc$OMP+        debut,lafin,rdayvrai,heure,dtphys,zplev,zplay,
530cc$OMP+        zphi,zphis,presnivs,clesphy0,zufi,zvfi,ztfi,
531cc$OMP+        zqfi,pvervel,zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
[630]532
[764]533c PRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
534c c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
535c c$OMP+                 zqfi_omp,pvervel_omp,zdufi_omp,zdvfi_omp,
536c c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp)
537
538c$OMP BARRIER
539      if (first_omp) then
[774]540        klon=klon_omp
[764]541
542        allocate(zplev_omp(klon,llm+1))
543        allocate(zplay_omp(klon,llm))
544        allocate(zphi_omp(klon,llm))
545        allocate(zphis_omp(klon))
546        allocate(presnivs_omp(llm))
547        allocate(zufi_omp(klon,llm))
548        allocate(zvfi_omp(klon,llm))
549        allocate(ztfi_omp(klon,llm))
[1114]550        allocate(zqfi_omp(klon,llm,nqtot))
[960]551c        allocate(pvervel_omp(klon,llm))
[764]552        allocate(zdufi_omp(klon,llm))
553        allocate(zdvfi_omp(klon,llm))
554        allocate(zdtfi_omp(klon,llm))
[1114]555        allocate(zdqfi_omp(klon,llm,nqtot))
[764]556        allocate(zdpsrf_omp(klon))
[985]557        allocate(flxwfi_omp(klon,llm))
[764]558        first_omp=.false.
559      endif
560       
561           
[774]562      klon=klon_omp
563      offset=klon_omp_begin-1
[764]564     
565      do l=1,llm+1
566        do i=1,klon
567          zplev_omp(i,l)=zplev(offset+i,l)
568        enddo
569      enddo
570         
571       do l=1,llm
572        do i=1,klon 
573          zplay_omp(i,l)=zplay(offset+i,l)
574        enddo
575      enddo
576       
577      do l=1,llm
578        do i=1,klon
579          zphi_omp(i,l)=zphi(offset+i,l)
580        enddo
581      enddo
582       
583      do i=1,klon
584        zphis_omp(i)=zphis(offset+i)
585      enddo
586     
587       
588      do l=1,llm
589        presnivs_omp(l)=presnivs(l)
590      enddo
591       
592      do l=1,llm
593        do i=1,klon
594          zufi_omp(i,l)=zufi(offset+i,l)
595        enddo
596      enddo
597       
598      do l=1,llm
599        do i=1,klon
600          zvfi_omp(i,l)=zvfi(offset+i,l)
601        enddo
602      enddo
603       
604      do l=1,llm
605        do i=1,klon
606          ztfi_omp(i,l)=ztfi(offset+i,l)
607        enddo
608      enddo
609       
[1114]610      do iq=1,nqtot
[764]611        do l=1,llm
612          do i=1,klon
613            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
614          enddo
615        enddo
616      enddo
617       
[960]618c      do l=1,llm
619c        do i=1,klon
620c         pvervel_omp(i,l)=pvervel(offset+i,l)
621c       enddo
622c      enddo
[764]623       
624      do l=1,llm
625        do i=1,klon
626          zdufi_omp(i,l)=zdufi(offset+i,l)
627        enddo
628      enddo
629       
630      do l=1,llm
631        do i=1,klon
632          zdvfi_omp(i,l)=zdvfi(offset+i,l)
633        enddo
634      enddo
635       
636      do l=1,llm
637        do i=1,klon
638          zdtfi_omp(i,l)=zdtfi(offset+i,l)
639        enddo
640      enddo
641       
[1114]642      do iq=1,nqtot
[764]643        do l=1,llm
644          do i=1,klon
645            zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
646          enddo
647        enddo
648      enddo
649       
650      do i=1,klon
651        zdpsrf_omp(i)=zdpsrf(offset+i)
652      enddo
[985]653
654      do l=1,llm
655        do i=1,klon
656          flxwfi_omp(i,l)=flxwfi(offset+i,l)
657        enddo
658      enddo
[764]659     
660c$OMP BARRIER
661cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
662     
[630]663      CALL physiq (klon,
664     .             llm,
665     .             debut,
666     .             lafin,
667     .             rdayvrai,
668     .             heure,
669     .             dtphys,
[764]670     .             zplev_omp,
671     .             zplay_omp,
672     .             zphi_omp,
673     .             zphis_omp,
674     .             presnivs_omp,
[630]675     .             clesphy0,
[764]676     .             zufi_omp,
677     .             zvfi_omp,
678     .             ztfi_omp,
679     .             zqfi_omp,
[960]680c     .             pvervel_omp,
681c#ifdef INCA
[985]682     .             flxwfi_omp,
[960]683c#endif
[764]684     .             zdufi_omp,
685     .             zdvfi_omp,
686     .             zdtfi_omp,
687     .             zdqfi_omp,
688     .             zdpsrf_omp,
689cIM diagnostique PVteta, Amip2         
690     .             pducov,
691     .             PVteta)
[630]692
[764]693cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
694
695c$OMP BARRIER
696
697      do l=1,llm+1
698        do i=1,klon
699          zplev(offset+i,l)=zplev_omp(i,l)
700        enddo
701      enddo
702         
703       do l=1,llm
704        do i=1,klon 
705          zplay(offset+i,l)=zplay_omp(i,l)
706        enddo
707      enddo
708       
709      do l=1,llm
710        do i=1,klon
711          zphi(offset+i,l)=zphi_omp(i,l)
712        enddo
713      enddo
714       
715
716      do i=1,klon
717        zphis(offset+i)=zphis_omp(i)
718      enddo
719     
720       
721      do l=1,llm
722        presnivs(l)=presnivs_omp(l)
723      enddo
724       
725      do l=1,llm
726        do i=1,klon
727          zufi(offset+i,l)=zufi_omp(i,l)
728        enddo
729      enddo
730       
731      do l=1,llm
732        do i=1,klon
733          zvfi(offset+i,l)=zvfi_omp(i,l)
734        enddo
735      enddo
736       
737      do l=1,llm
738        do i=1,klon
739          ztfi(offset+i,l)=ztfi_omp(i,l)
740        enddo
741      enddo
742       
[1114]743      do iq=1,nqtot
[764]744        do l=1,llm
745          do i=1,klon
746            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
747          enddo
748        enddo
749      enddo
750       
[960]751c      do l=1,llm
752c        do i=1,klon
753c         pvervel(offset+i,l)=pvervel_omp(i,l)
754c       enddo
755c      enddo
[764]756       
757      do l=1,llm
758        do i=1,klon
759          zdufi(offset+i,l)=zdufi_omp(i,l)
760        enddo
761      enddo
762       
763      do l=1,llm
764        do i=1,klon
765          zdvfi(offset+i,l)=zdvfi_omp(i,l)
766        enddo
767      enddo
768       
769      do l=1,llm
770        do i=1,klon
771          zdtfi(offset+i,l)=zdtfi_omp(i,l)
772        enddo
773      enddo
774       
[1114]775      do iq=1,nqtot
[764]776        do l=1,llm
777          do i=1,klon
778            zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
779          enddo
780        enddo
781      enddo
782       
783      do i=1,klon
784        zdpsrf(offset+i)=zdpsrf_omp(i)
785      enddo
786     
787
788cc$OMP END PARALLEL
789      klon=klon_mpi
[630]790500   CONTINUE
[764]791c$OMP BARRIER
[630]792
[764]793c$OMP MASTER
794cym      call WriteField_phy('zdtfi',zdtfi(:,:),llm)
[630]795      call stop_timer(timer_physic)
[764]796c$OMP END MASTER
[1000]797
798      IF (using_mpi) THEN
799           
[630]800      if (MPI_rank>0) then
[764]801
802c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
803       DO l=1,llm     
804        du_send(1:iim,l)=zdufi(1:iim,l)
805        dv_send(1:iim,l)=zdvfi(1:iim,l)
806       ENDDO
807c$OMP END DO NOWAIT       
808
809c$OMP BARRIER
[1000]810#ifdef CPP_MPI
[764]811c$OMP MASTER
[985]812!$OMP CRITICAL (MPI)
[630]813        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
[764]814     &                   COMM_LMDZ,Req(1),ierr)
[630]815        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
[764]816     &                  COMM_LMDZ,Req(2),ierr)
[985]817!$OMP END CRITICAL (MPI)
[764]818c$OMP END MASTER
[1000]819#endif
[764]820c$OMP BARRIER
[630]821     
822      endif
823   
824      if (MPI_rank<MPI_Size-1) then
[764]825c$OMP BARRIER
[1000]826#ifdef CPP_MPI
[764]827c$OMP MASTER     
[985]828!$OMP CRITICAL (MPI)
[630]829        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
[764]830     &                 COMM_LMDZ,Req(3),ierr)
[630]831        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
[764]832     &                 COMM_LMDZ,Req(4),ierr)
[985]833!$OMP END CRITICAL (MPI)
[764]834c$OMP END MASTER
[1000]835#endif
[630]836      endif
[764]837
838c$OMP BARRIER
[1000]839
840
841#ifdef CPP_MPI
[764]842c$OMP MASTER   
[985]843!$OMP CRITICAL (MPI)
[630]844      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
845        call MPI_WAITALL(4,Req(1),Status,ierr)
846      else if (MPI_rank>0) then
847        call MPI_WAITALL(2,Req(1),Status,ierr)
848      else if (MPI_rank <MPI_Size-1) then
849        call MPI_WAITALL(2,Req(3),Status,ierr)
850      endif
[985]851!$OMP END CRITICAL (MPI)
[764]852c$OMP END MASTER
[1000]853#endif
854
[764]855c$OMP BARRIER     
856
[1000]857      ENDIF ! using_mpi
858     
859     
[764]860c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
861      DO l=1,llm
862           
863        zdufi2(1:klon,l)=zdufi(1:klon,l)
864        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
865           
866        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
867        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
868 
[774]869         pdhfi(:,jj_begin,l)=0
870         pdqfi(:,jj_begin,l,:)=0
871         pdufi(:,jj_begin,l)=0
872         pdvfi(:,jj_begin,l)=0
[764]873         
[774]874         if (.not. is_south_pole) then
875           pdhfi(:,jj_end,l)=0
876           pdqfi(:,jj_end,l,:)=0
877           pdufi(:,jj_end,l)=0
878           pdvfi(:,jj_end,l)=0
[764]879         endif
[630]880     
[764]881       ENDDO
882c$OMP END DO NOWAIT
[630]883
[764]884c$OMP MASTER
[774]885       pdpsfi(:,jj_begin)=0   
886       if (.not. is_south_pole) then
887         pdpsfi(:,jj_end)=0
[630]888       endif
[764]889c$OMP END MASTER
[630]890c-----------------------------------------------------------------------
891c   transformation des tendances physiques en tendances dynamiques:
892c   ---------------------------------------------------------------
893
894c  tendance sur la pression :
895c  -----------------------------------
896      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
897c
898c   62. enthalpie potentielle
899c   ---------------------
900     
901      kstart=1
902      kend=klon
903
[774]904      if (is_north_pole) kstart=2
905      if (is_south_pole)  kend=klon-1
[630]906
[764]907c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]908      DO l=1,llm
909
[764]910!!cdir NODEP
[630]911        do ig0=kstart,kend
[774]912          i=index_i(ig0)
913          j=index_j(ig0)
[630]914          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
915          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
916         enddo         
917
[774]918        if (is_north_pole) then
[630]919            DO i=1,iip1
920              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
921            enddo
922        endif
923       
[774]924        if (is_south_pole) then
[630]925            DO i=1,iip1
926              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
927            ENDDO
928        endif
929      ENDDO
[764]930c$OMP END DO NOWAIT
[630]931     
932c   62. humidite specifique
933c   ---------------------
934
[1114]935      DO iq=1,nqtot
[764]936c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]937         DO l=1,llm
[764]938!!cdir NODEP
[630]939           do ig0=kstart,kend
[774]940             i=index_i(ig0)
941             j=index_j(ig0)
[630]942             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
943             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
944           enddo
945           
[774]946           if (is_north_pole) then
[630]947             do i=1,iip1
948               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
949             enddo
950           endif
951           
[774]952           if (is_south_pole) then
[630]953             do i=1,iip1
954               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
955             enddo
956           endif
957           
958         ENDDO
[764]959c$OMP END DO NOWAIT
[630]960      ENDDO
961
962c   63. traceurs
963c   ------------
964C     initialisation des tendances
[764]965
966c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
967      DO l=1,llm
968        pdqfi(:,:,l,:)=0.
969      ENDDO
970c$OMP END DO NOWAIT     
971
[630]972C
973
[1114]974      DO iq=1,nqtot
[630]975         iiq=niadv(iq)
[764]976c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]977         DO l=1,llm
978
[764]979!!cdir NODEP           
[630]980             DO ig0=kstart,kend
[774]981              i=index_i(ig0)
982              j=index_j(ig0)
[630]983              pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
984              if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
985            ENDDO
986           
[774]987            IF (is_north_pole) then
[630]988              DO i=1,iip1
989                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
990              ENDDO
991            ENDIF
992           
[774]993            IF (is_south_pole) then
[630]994              DO i=1,iip1
995                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
996              ENDDO
997            ENDIF
998           
999         ENDDO
[764]1000c$OMP END DO NOWAIT     
[630]1001      ENDDO
1002     
1003c   65. champ u:
1004c   ------------
[764]1005c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]1006      DO l=1,llm
[764]1007!!cdir NODEP
[630]1008         do ig0=kstart,kend
[774]1009           i=index_i(ig0)
1010           j=index_j(ig0)
[630]1011           
1012           if (i/=iim) then
1013             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1014           endif
1015           
1016           if (i==1) then
1017              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
1018     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
[764]1019             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
[630]1020           endif
1021         
1022         enddo
1023         
[774]1024         if (is_north_pole) then
[630]1025           DO i=1,iip1
1026            pdufi(i,1,l)    = 0.
1027           ENDDO
1028         endif
1029         
[774]1030         if (is_south_pole) then
[630]1031           DO i=1,iip1
1032            pdufi(i,jjp1,l) = 0.
1033           ENDDO
1034         endif
1035         
1036      ENDDO
[764]1037c$OMP END DO NOWAIT
[630]1038
1039c   67. champ v:
1040c   ------------
1041
1042      kstart=1
1043      kend=klon
1044
[774]1045      if (is_north_pole) kstart=2
1046      if (is_south_pole)  kend=klon-1-iim
[630]1047     
[764]1048c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]1049      DO l=1,llm
[764]1050!!cdir NODEP
[630]1051        do ig0=kstart,kend
[774]1052           i=index_i(ig0)
1053           j=index_j(ig0)
[630]1054           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
1055           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1056     $                                      zdvfi2(ig0+iim,l))
1057     $                                    *cv(i,j)
1058        enddo
1059         
1060      ENDDO
[764]1061c$OMP END DO NOWAIT
[630]1062
1063
1064c   68. champ v pres des poles:
1065c   ---------------------------
1066c      v = U * cos(long) + V * SIN(long)
1067
[774]1068      if (is_north_pole) then
[764]1069
1070c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]1071        DO l=1,llm
1072
1073          DO i=1,iim
1074            pdvfi(i,1,l)=
1075     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
1076       
1077            pdvfi(i,1,l)=
1078     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
1079          ENDDO
1080
1081          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
1082
1083        ENDDO
[764]1084c$OMP END DO NOWAIT
[630]1085
1086      endif   
1087     
[774]1088      if (is_south_pole) then
[764]1089
1090c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1091         DO l=1,llm
[630]1092 
1093           DO i=1,iim
1094              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1095     $        +zdvfi(klon,l)*SIN(rlonv(i))
1096
1097              pdvfi(i,jjm,l)=
1098     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
1099           ENDDO
1100
1101           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1102
1103        ENDDO
[764]1104c$OMP END DO NOWAIT
[630]1105     
1106      endif
1107c-----------------------------------------------------------------------
1108
1109700   CONTINUE
1110 
1111      firstcal = .FALSE.
1112
1113      RETURN
1114      END
Note: See TracBrowser for help on using the repository browser.