source: LMDZ4/branches/V3_test/libf/dyn3dpar/calfis_p.F @ 5071

Last change on this file since 5071 was 709, checked in by Laurent Fairhead, 18 years ago

Nouvelles versions de la dynamique YM
LF

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