source: LMDZ5/trunk/libf/calfis_loc.F @ 1630

Last change on this file since 1630 was 1630, checked in by Laurent Fairhead, 12 years ago

Importation initiale du répertoire dyn3dmem

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