source: trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F @ 1022

Last change on this file since 1022 was 1019, checked in by emillour, 13 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

  • General stuff:
  • makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

  • bld.cfg: adaptations to enable compiling using multiple threads
  • build_gcm: adaptations to enable compiling using multiple threads
  • makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change
  • bibio:
  • wxios.F90 : Added for possible future use of XIOS library
  • filtrez:
  • mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
  • filtreg_mod : limit use of FFT to parallel mode
  • mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz
  • dyn3d:
  • abort_gcm.F : add things for xios
  • advtrac.F90 : minor change in CFL outputs
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • inigeom.F : move two computations outside loop
  • dyn3dpar:
  • parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • leapfrog_p.F : add xios + correction for times in Newtonian case
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • inigeom.F : move two computations outside loop

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