source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/calfis_p.F @ 2407

Last change on this file since 2407 was 1407, checked in by Laurent Fairhead, 14 years ago

Some extraneous prints deleted for optimisation purposes

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