source: LMDZ4/tags/LMDZ4_V3_3/libf/dyn3dpar/calfis_p.F @ 1401

Last change on this file since 1401 was 774, checked in by Laurent Fairhead, 17 years ago

Suite du merge entre la version et la HEAD: quelques modifications de
Yann sur le

LF

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