source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/dynphy_lonlat/calfis_p.F @ 5409

Last change on this file since 5409 was 2604, checked in by Ehouarn Millour, 8 years ago

Add Exner function to the call_physiq arguments (not used by the Earth physics) to harmonize physics/dynamics interface.
EM

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