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

Last change on this file since 1142 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
Line 
1!
2! $Header$
3!
4C
5C
6      SUBROUTINE calfis_p(lafin,
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
34      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
35      USE parallel, ONLY : omp_chunk, using_mpi
36      USE mod_interface_dyn_phys
37      USE Write_Field
38      Use Write_field_p
39      USE Times
40      USE IOPHY
41      USE infotrac
42
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
102      INTEGER ngridmx
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"
109#ifdef CPP_MPI
110      include 'mpif.h'
111#endif
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)
121      REAL pq(iip1,jjp1,llm,nqtot)
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)
128      REAL pdq(iip1,jjp1,llm,nqtot)
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)
137      REAL pdqfi(iip1,jjp1,llm,nqtot)
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
149      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
150      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
151      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
152c
153      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
154      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
155c
156      REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
157      REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
158c
159c      REAL,ALLOCATABLE,SAVE :: pvervel(:,:)
160c
161      REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
162      REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
163      REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
164      REAL,SAVE,ALLOCATABLE ::  flxwfi(:,:)     ! Flux de masse verticale sur la grille physiq
165
166c
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(:,:,:)
176c      REAL,ALLOCATABLE,SAVE :: pvervel_omp(:,:)
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(:)
182      REAL,SAVE,ALLOCATABLE ::  flxwfi_omp(:,:)     ! Flux de masse verticale sur la grille physiq
183
184c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
185c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
186c$OMP+                 zqfi_omp,zdufi_omp,zdvfi_omp,
187c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp)       
188
189      LOGICAL,SAVE :: first_omp=.true.
190c$OMP THREADPRIVATE(first_omp)
191     
192      REAL zsin(iim),zcos(iim),z1(iim)
193      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
194      REAL unskap, pksurcp
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     
203      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
204     
205      REAL SSUM
206
207      LOGICAL firstcal, debut
208      DATA firstcal/.true./
209      SAVE firstcal,debut
210c$OMP THREADPRIVATE(firstcal,debut)
211      REAL rdayvrai
212     
213      REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
214      INTEGER :: ierr
215#ifdef CPP_MPI
216      INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
217#else
218      INTEGER,dimension(1,4) :: Status
219#endif
220      INTEGER, dimension(4) :: Req
221      REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
222      integer :: k,kstart,kend
223      INTEGER :: offset 
224c
225c-----------------------------------------------------------------------
226c
227c    1. Initialisations :
228c    --------------------
229c
230
231      klon=klon_mpi
232     
233      PVteta(:,:)=0.
234           
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.
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))
255      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
256      ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
257      ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
258c      ALLOCATE(pvervel(klon,llm))
259      ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
260      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
261      ALLOCATE(zdpsrf(klon))
262      ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
263      ALLOCATE(flxwfi(klon,llm))
264c$OMP END MASTER
265c$OMP BARRIER     
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
279c$OMP MASTER
280      call start_timer(timer_physic)
281c$OMP END MASTER
282
283c$OMP MASTER             
284      do ig0=1,klon
285        i=index_i(ig0)
286        j=index_j(ig0)
287        zpsrf(ig0)=pps(i,j)
288      enddo
289c$OMP END MASTER
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
303c      print *,omp_rank,'klon--->',klon
304c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
305      DO l = 1, llmp1
306        do ig0=1,klon
307          i=index_i(ig0)
308          j=index_j(ig0)
309          zplev( ig0,l ) = pp(i,j,l)
310        enddo
311      ENDDO
312c$OMP END DO NOWAIT
313c
314c
315
316c   43. temperature naturelle (en K) et pressions milieux couches .
317c   ---------------------------------------------------------------
318c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
319      DO l=1,llm
320
321        do ig0=1,klon
322          i=index_i(ig0)
323          j=index_j(ig0)
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
331c$OMP END DO NOWAIT
332
333c   43.bis traceurs
334c   ---------------
335c
336
337      DO iq=1,nqtot
338         iiq=niadv(iq)
339c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
340         DO l=1,llm
341           do ig0=1,klon
342             i=index_i(ig0)
343             j=index_j(ig0)
344             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
345           enddo
346         ENDDO
347c$OMP END DO NOWAIT     
348      ENDDO
349
350c   convergence dynamique pour les traceurs "EAU"
351
352      DO iq=1,2
353c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
354         DO l=1,llm
355           do ig0=1,klon
356             i=index_i(ig0)
357             j=index_j(ig0)
358c             pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
359           enddo
360         ENDDO
361c$OMP END DO NOWAIT     
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)
370
371      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
372
373c$OMP BARRIER
374
375c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
376      DO l=1,llm
377         DO ig=1,klon
378           zphi(ig,l)=zphi(ig,l)-zphis(ig)
379         ENDDO
380      ENDDO
381c$OMP END DO NOWAIT
382     
383c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
384c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
385c      de masse est calclue dans advtrac_p.F 
386c
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
398
399c
400c   45. champ u:
401c   ------------
402
403      kstart=1
404      kend=klon
405     
406      if (is_north_pole) kstart=2
407      if (is_south_pole) kend=klon-1
408     
409c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
410      DO l=1,llm
411        do ig0=kstart,kend
412          i=index_i(ig0)
413          j=index_j(ig0)
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
427c$OMP END DO NOWAIT
428c   46.champ v:
429c   -----------
430c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
431      DO l=1,llm
432        DO ig0=kstart,kend
433          i=index_i(ig0)
434          j=index_j(ig0)
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
442c$OMP END DO NOWAIT
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
449      if (is_north_pole) then
450c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
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
473c$OMP END DO NOWAIT     
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
482      if (is_south_pole) then
483c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
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
506c$OMP END DO NOWAIT       
507      endif
508
509
510      IF (is_sequential) THEN
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
518
519c On change de grille, dynamique vers physiq, pour le flux de masse verticale
520      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
521
522c-----------------------------------------------------------------------
523c   Appel de la physique:
524c   ---------------------
525
526cc$OMP  PARALLEL DEFAULT(NONE)
527cc$OMP+ PRIVATE(i,l,offset,iq)
528cc$OMP+ SHARED(klon_omp_nb,nqtot,klon_omp_begin,
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)
532
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
540        klon=klon_omp
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))
550        allocate(zqfi_omp(klon,llm,nqtot))
551c        allocate(pvervel_omp(klon,llm))
552        allocate(zdufi_omp(klon,llm))
553        allocate(zdvfi_omp(klon,llm))
554        allocate(zdtfi_omp(klon,llm))
555        allocate(zdqfi_omp(klon,llm,nqtot))
556        allocate(zdpsrf_omp(klon))
557        allocate(flxwfi_omp(klon,llm))
558        first_omp=.false.
559      endif
560       
561           
562      klon=klon_omp
563      offset=klon_omp_begin-1
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       
610      do iq=1,nqtot
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       
618c      do l=1,llm
619c        do i=1,klon
620c         pvervel_omp(i,l)=pvervel(offset+i,l)
621c       enddo
622c      enddo
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       
642      do iq=1,nqtot
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
653
654      do l=1,llm
655        do i=1,klon
656          flxwfi_omp(i,l)=flxwfi(offset+i,l)
657        enddo
658      enddo
659     
660c$OMP BARRIER
661cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
662     
663      CALL physiq (klon,
664     .             llm,
665     .             debut,
666     .             lafin,
667     .             rdayvrai,
668     .             heure,
669     .             dtphys,
670     .             zplev_omp,
671     .             zplay_omp,
672     .             zphi_omp,
673     .             zphis_omp,
674     .             presnivs_omp,
675     .             clesphy0,
676     .             zufi_omp,
677     .             zvfi_omp,
678     .             ztfi_omp,
679     .             zqfi_omp,
680c     .             pvervel_omp,
681c#ifdef INCA
682     .             flxwfi_omp,
683c#endif
684     .             zdufi_omp,
685     .             zdvfi_omp,
686     .             zdtfi_omp,
687     .             zdqfi_omp,
688     .             zdpsrf_omp,
689cIM diagnostique PVteta, Amip2         
690     .             pducov,
691     .             PVteta)
692
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       
743      do iq=1,nqtot
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       
751c      do l=1,llm
752c        do i=1,klon
753c         pvervel(offset+i,l)=pvervel_omp(i,l)
754c       enddo
755c      enddo
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       
775      do iq=1,nqtot
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
790500   CONTINUE
791c$OMP BARRIER
792
793c$OMP MASTER
794cym      call WriteField_phy('zdtfi',zdtfi(:,:),llm)
795      call stop_timer(timer_physic)
796c$OMP END MASTER
797
798      IF (using_mpi) THEN
799           
800      if (MPI_rank>0) then
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
810#ifdef CPP_MPI
811c$OMP MASTER
812!$OMP CRITICAL (MPI)
813        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
814     &                   COMM_LMDZ,Req(1),ierr)
815        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
816     &                  COMM_LMDZ,Req(2),ierr)
817!$OMP END CRITICAL (MPI)
818c$OMP END MASTER
819#endif
820c$OMP BARRIER
821     
822      endif
823   
824      if (MPI_rank<MPI_Size-1) then
825c$OMP BARRIER
826#ifdef CPP_MPI
827c$OMP MASTER     
828!$OMP CRITICAL (MPI)
829        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
830     &                 COMM_LMDZ,Req(3),ierr)
831        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
832     &                 COMM_LMDZ,Req(4),ierr)
833!$OMP END CRITICAL (MPI)
834c$OMP END MASTER
835#endif
836      endif
837
838c$OMP BARRIER
839
840
841#ifdef CPP_MPI
842c$OMP MASTER   
843!$OMP CRITICAL (MPI)
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
851!$OMP END CRITICAL (MPI)
852c$OMP END MASTER
853#endif
854
855c$OMP BARRIER     
856
857      ENDIF ! using_mpi
858     
859     
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 
869         pdhfi(:,jj_begin,l)=0
870         pdqfi(:,jj_begin,l,:)=0
871         pdufi(:,jj_begin,l)=0
872         pdvfi(:,jj_begin,l)=0
873         
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
879         endif
880     
881       ENDDO
882c$OMP END DO NOWAIT
883
884c$OMP MASTER
885       pdpsfi(:,jj_begin)=0   
886       if (.not. is_south_pole) then
887         pdpsfi(:,jj_end)=0
888       endif
889c$OMP END MASTER
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
904      if (is_north_pole) kstart=2
905      if (is_south_pole)  kend=klon-1
906
907c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
908      DO l=1,llm
909
910!!cdir NODEP
911        do ig0=kstart,kend
912          i=index_i(ig0)
913          j=index_j(ig0)
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
918        if (is_north_pole) then
919            DO i=1,iip1
920              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
921            enddo
922        endif
923       
924        if (is_south_pole) then
925            DO i=1,iip1
926              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
927            ENDDO
928        endif
929      ENDDO
930c$OMP END DO NOWAIT
931     
932c   62. humidite specifique
933c   ---------------------
934
935      DO iq=1,nqtot
936c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
937         DO l=1,llm
938!!cdir NODEP
939           do ig0=kstart,kend
940             i=index_i(ig0)
941             j=index_j(ig0)
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           
946           if (is_north_pole) then
947             do i=1,iip1
948               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
949             enddo
950           endif
951           
952           if (is_south_pole) then
953             do i=1,iip1
954               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
955             enddo
956           endif
957           
958         ENDDO
959c$OMP END DO NOWAIT
960      ENDDO
961
962c   63. traceurs
963c   ------------
964C     initialisation des tendances
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
972C
973
974      DO iq=1,nqtot
975         iiq=niadv(iq)
976c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
977         DO l=1,llm
978
979!!cdir NODEP           
980             DO ig0=kstart,kend
981              i=index_i(ig0)
982              j=index_j(ig0)
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           
987            IF (is_north_pole) then
988              DO i=1,iip1
989                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
990              ENDDO
991            ENDIF
992           
993            IF (is_south_pole) then
994              DO i=1,iip1
995                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
996              ENDDO
997            ENDIF
998           
999         ENDDO
1000c$OMP END DO NOWAIT     
1001      ENDDO
1002     
1003c   65. champ u:
1004c   ------------
1005c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1006      DO l=1,llm
1007!!cdir NODEP
1008         do ig0=kstart,kend
1009           i=index_i(ig0)
1010           j=index_j(ig0)
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)
1019             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1020           endif
1021         
1022         enddo
1023         
1024         if (is_north_pole) then
1025           DO i=1,iip1
1026            pdufi(i,1,l)    = 0.
1027           ENDDO
1028         endif
1029         
1030         if (is_south_pole) then
1031           DO i=1,iip1
1032            pdufi(i,jjp1,l) = 0.
1033           ENDDO
1034         endif
1035         
1036      ENDDO
1037c$OMP END DO NOWAIT
1038
1039c   67. champ v:
1040c   ------------
1041
1042      kstart=1
1043      kend=klon
1044
1045      if (is_north_pole) kstart=2
1046      if (is_south_pole)  kend=klon-1-iim
1047     
1048c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1049      DO l=1,llm
1050!!cdir NODEP
1051        do ig0=kstart,kend
1052           i=index_i(ig0)
1053           j=index_j(ig0)
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
1061c$OMP END DO NOWAIT
1062
1063
1064c   68. champ v pres des poles:
1065c   ---------------------------
1066c      v = U * cos(long) + V * SIN(long)
1067
1068      if (is_north_pole) then
1069
1070c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
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
1084c$OMP END DO NOWAIT
1085
1086      endif   
1087     
1088      if (is_south_pole) then
1089
1090c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1091         DO l=1,llm
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
1104c$OMP END DO NOWAIT
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.