source: LMDZ5/trunk/libf/dynphy_lonlat/calfis_p.F @ 2526

Last change on this file since 2526 was 2429, checked in by Laurent Fairhead, 9 years ago

Correction on the calculation of the surface of the grid at the poles (problem was introduced
in r2222).
Due to the different distribution of OMP tasks in the dynamics and the physics, had to
introduce 2 new logical variables, is_pole_north_phy and is_pole_south_phy, and so decided
to rename the old is_north_pole/is_south_pole to is_north_pole_dyn/is_south_pole_dyn to
stay coherent and, hopefully, clear things up a bit.
LF

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