source: LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/calfis_loc.F @ 3400

Last change on this file since 3400 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:keywords set to Id
File size: 32.8 KB
Line 
1!
2! $Id: calfis_loc.F 2604 2016-07-26 15:37:18Z fairhead $
3!
4C
5C
6      SUBROUTINE calfis_loc(lafin,
7     $                  jD_cur, jH_cur,
8     $                  pucov,
9     $                  pvcov,
10     $                  pteta,
11     $                  pq,
12     $                  pmasse,
13     $                  pps,
14     $                  pp,
15     $                  ppk,
16     $                  pphis,
17     $                  pphi,
18     $                  pducov,
19     $                  pdvcov,
20     $                  pdteta,
21     $                  pdq,
22     $                  flxw,
23     $                  pdufi,
24     $                  pdvfi,
25     $                  pdhfi,
26     $                  pdqfi,
27     $                  pdpsfi)
28#ifdef CPP_PHYS
29! If using physics
30c
31c    Auteur :  P. Le Van, F. Hourdin
32c   .........
33      USE dimphy
34      USE mod_phys_lmdz_mpi_data, mpi_root_xx=>mpi_master
35      USE mod_phys_lmdz_omp_data, ONLY: klon_omp, klon_omp_begin
36      USE mod_const_mpi, ONLY: COMM_LMDZ
37      USE mod_interface_dyn_phys
38      USE IOPHY
39#endif
40#ifdef CPP_PARA
41      USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
42     $                        ,jj_begin_dyn=>jj_begin,jj_end_dyn=>jj_end
43      USE Write_Field
44      Use Write_field_p
45      USE Times
46#endif
47      USE infotrac, ONLY: nqtot, niadv, tname
48      USE control_mod, ONLY: planet_type, nsplit_phys
49#ifdef CPP_PHYS
50      USE callphysiq_mod, ONLY: call_physiq
51#endif
52      USE comvert_mod, ONLY: preff, presnivs
53      USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
54
55#ifdef CPP_PARA
56      IMPLICIT NONE
57c=======================================================================
58c
59c   1. rearrangement des tableaux et transformation
60c      variables dynamiques  >  variables physiques
61c   2. calcul des termes physiques
62c   3. retransformation des tendances physiques en tendances dynamiques
63c
64c   remarques:
65c   ----------
66c
67c    - les vents sont donnes dans la physique par leurs composantes
68c      naturelles.
69c    - la variable thermodynamique de la physique est une variable
70c      intensive :   T
71c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
72c    - les deux seules variables dependant de la geometrie necessaires
73c      pour la physique sont la latitude pour le rayonnement et
74c      l'aire de la maille quand on veut integrer une grandeur
75c      horizontalement.
76c    - les points de la physique sont les points scalaires de la
77c      la dynamique; numerotation:
78c          1 pour le pole nord
79c          (jjm-1)*iim pour l'interieur du domaine
80c          ngridmx pour le pole sud
81c      ---> ngridmx=2+(jjm-1)*iim
82c
83c     Input :
84c     -------
85c       ecritphy        frequence d'ecriture (en jours)de histphy
86c       pucov           covariant zonal velocity
87c       pvcov           covariant meridional velocity
88c       pteta           potential temperature
89c       pps             surface pressure
90c       pmasse          masse d'air dans chaque maille
91c       pts             surface temperature  (K)
92c       callrad         clef d'appel au rayonnement
93c
94c    Output :
95c    --------
96c        pdufi          tendency for the natural zonal velocity (ms-1)
97c        pdvfi          tendency for the natural meridional velocity
98c        pdhfi          tendency for the potential temperature
99c        pdtsfi         tendency for the surface temperature
100c
101c        pdtrad         radiative tendencies  \  both input
102c        pfluxrad       radiative fluxes      /  and output
103c
104c=======================================================================
105c
106c-----------------------------------------------------------------------
107c
108c    0.  Declarations :
109c    ------------------
110
111      include "dimensions.h"
112      include "paramet.h"
113
114      INTEGER ngridmx
115      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
116
117      include "comgeom2.h"
118      include "iniprint.h"
119#ifdef CPP_MPI
120      include 'mpif.h'
121#endif
122c    Arguments :
123c    -----------
124      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
125      REAL,INTENT(IN):: jD_cur, jH_cur
126      REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
127      REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
128      REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
129      REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
130      REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
131      REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
132      REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential
133
134      REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
135      REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
136      REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
137      REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used
138
139      REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
140      REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
141      REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
142      REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
143
144      ! tendencies (in */s) from the physics
145      REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
146      REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
147      REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
148      REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
149      REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)
150
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(:,:,:)
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
387c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
388         DO l=1,llm
389!CDIR ON_ADB(index_i)
390!CDIR ON_ADB(index_j)
391           do ig0=1,klon
392             i=index_i(ig0)
393             j=index_j(ig0)
394             zphi(ig0,l)  = pphi(i,j,l)
395           enddo
396         ENDDO
397c$OMP END DO NOWAIT         
398
399c      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
400
401c$OMP MASTER
402!CDIR ON_ADB(index_i)
403!CDIR ON_ADB(index_j)
404           do ig0=1,klon
405             i=index_i(ig0)
406             j=index_j(ig0)
407             zphis(ig0)  = pphis(i,j)
408           enddo
409c$OMP END MASTER
410
411
412c      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
413
414c$OMP BARRIER
415
416c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
417      DO l=1,llm
418         DO ig=1,klon
419           zphi(ig,l)=zphi(ig,l)-zphis(ig)
420         ENDDO
421      ENDDO
422c$OMP END DO NOWAIT
423     
424
425c
426c   45. champ u:
427c   ------------
428
429      kstart=1
430      kend=klon
431     
432      if (is_north_pole_dyn) kstart=2
433      if (is_south_pole_dyn) kend=klon-1
434     
435c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
436      DO l=1,llm
437!CDIR ON_ADB(index_i)
438!CDIR ON_ADB(index_j)
439!CDIR SPARSE
440        do ig0=kstart,kend
441          i=index_i(ig0)
442          j=index_j(ig0)
443          if (i==1) then
444            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
445     $                         + pucov(1,j,l)/cu(1,j) )
446          else
447            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j)
448     $                       + pucov(i,j,l)/cu(i,j) )
449          endif
450        enddo
451      ENDDO
452c$OMP END DO NOWAIT
453
454c
455C  Alvaro de la Camara (May 2014)
456C  46.1 Calcul de la vorticite et passage sur la grille physique
457C  --------------------------------------------------------------
458
459      jjb=jj_begin_dyn-1
460      jje=jj_end_dyn+1
461      if (is_north_pole_dyn) jjb=1
462      if (is_south_pole_dyn) jje=jjm
463
464c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
465
466      DO l=1,llm
467        do i=1,iim
468          do j=jjb,jje
469            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
470     $                   + pucov(i,j+1,l) - pucov(i,j,l))
471     $                   / (cu(i,j)+cu(i,j+1))
472     $                   / (cv(i+1,j)+cv(i,j)) *4
473          enddo
474        enddo
475      ENDDO
476
477
478c   46.2champ v:
479c   -----------
480
481c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
482      DO l=1,llm
483!CDIR ON_ADB(index_i)
484!CDIR ON_ADB(index_j)
485        DO ig0=kstart,kend
486          i=index_i(ig0)
487          j=index_j(ig0)
488          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1)
489     $                       + pvcov(i,j,l)/cv(i,j) )
490          if (j==1 .OR. j==jjp1) then !  AdlC MAY 2014
491            zrfi(ig0,l) = 0 !  AdlC MAY 2014
492          else
493            if(i==1)then
494            zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
495     $                   +zrot(1,j-1,l)+zrot(1,j,l))   !  AdlC MAY 2014
496            else
497            zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
498     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
499            endif
500          endif
501
502   
503         ENDDO
504      ENDDO
505c$OMP END DO NOWAIT
506
507c   47. champs de vents aux pole nord   
508c   ------------------------------
509c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
510c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
511
512      if (is_north_pole_dyn) then
513c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
514        DO l=1,llm
515
516           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
517           DO i=2,iim
518              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
519           ENDDO
520 
521           DO i=1,iim
522              zcos(i)   = COS(rlonv(i))*z1(i)
523              zsin(i)   = SIN(rlonv(i))*z1(i)
524           ENDDO
525 
526           zufi(1,l)  = SSUM(iim,zcos,1)/pi
527           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
528           zrfi(1,l)  = 0.
529 
530        ENDDO
531c$OMP END DO NOWAIT     
532      endif
533
534
535c   48. champs de vents aux pole sud:
536c   ---------------------------------
537c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
538c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
539
540      if (is_south_pole_dyn) then
541c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
542        DO l=1,llm
543 
544         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
545           DO i=2,iim
546             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
547           ENDDO
548 
549           DO i=1,iim
550              zcos(i)    = COS(rlonv(i))*z1(i)
551              zsin(i)    = SIN(rlonv(i))*z1(i)
552           ENDDO
553 
554           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
555           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
556           zrfi(klon,l)  = 0.
557        ENDDO
558c$OMP END DO NOWAIT       
559      endif
560
561c On change de grille, dynamique vers physiq, pour le flux de masse verticale
562c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
563         DO l=1,llm
564!CDIR ON_ADB(index_i)
565!CDIR ON_ADB(index_j)
566           do ig0=1,klon
567             i=index_i(ig0)
568             j=index_j(ig0)
569             flxwfi(ig0,l)  = flxw(i,j,l)
570           enddo
571         ENDDO
572c$OMP END DO NOWAIT
573
574c      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
575
576c-----------------------------------------------------------------------
577c   Appel de la physique:
578c   ---------------------
579
580
581c$OMP BARRIER
582      if (first_omp) then
583        klon=klon_omp
584
585        allocate(zplev_omp(klon,llm+1))
586        allocate(zplay_omp(klon,llm))
587        allocate(zpk_omp(klon,llm))
588        allocate(zphi_omp(klon,llm))
589        allocate(zphis_omp(klon))
590        allocate(presnivs_omp(llm))
591        allocate(zufi_omp(klon,llm))
592        allocate(zvfi_omp(klon,llm))
593        allocate(zrfi_omp(klon,llm))  ! LG Ari 2014
594        allocate(ztfi_omp(klon,llm))
595        allocate(zqfi_omp(klon,llm,nqtot))
596        allocate(zdufi_omp(klon,llm))
597        allocate(zdvfi_omp(klon,llm))
598        allocate(zdtfi_omp(klon,llm))
599        allocate(zdqfi_omp(klon,llm,nqtot))
600        allocate(zdufic_omp(klon,llm))
601        allocate(zdvfic_omp(klon,llm))
602        allocate(zdtfic_omp(klon,llm))
603        allocate(zdqfic_omp(klon,llm,nqtot))
604        allocate(zdpsrf_omp(klon))
605        allocate(flxwfi_omp(klon,llm))
606        first_omp=.false.
607      endif
608       
609           
610      klon=klon_omp
611      offset=klon_omp_begin-1
612     
613      do l=1,llm+1
614        do i=1,klon
615          zplev_omp(i,l)=zplev(offset+i,l)
616        enddo
617      enddo
618         
619       do l=1,llm
620        do i=1,klon 
621          zplay_omp(i,l)=zplay(offset+i,l)
622        enddo
623      enddo
624       
625       do l=1,llm
626        do i=1,klon 
627          zpk_omp(i,l)=zpk(offset+i,l)
628        enddo
629      enddo
630       
631      do l=1,llm
632        do i=1,klon
633          zphi_omp(i,l)=zphi(offset+i,l)
634        enddo
635      enddo
636       
637      do i=1,klon
638        zphis_omp(i)=zphis(offset+i)
639      enddo
640     
641       
642      do l=1,llm
643        presnivs_omp(l)=presnivs(l)
644      enddo
645       
646      do l=1,llm
647        do i=1,klon
648          zufi_omp(i,l)=zufi(offset+i,l)
649        enddo
650      enddo
651       
652      do l=1,llm
653        do i=1,klon
654          zvfi_omp(i,l)=zvfi(offset+i,l)
655        enddo
656      enddo
657       
658      do l=1,llm
659        do i=1,klon
660          zrfi_omp(i,l)=zrfi(offset+i,l)
661        enddo
662      enddo
663       
664      do l=1,llm
665        do i=1,klon
666          ztfi_omp(i,l)=ztfi(offset+i,l)
667        enddo
668      enddo
669       
670      do iq=1,nqtot
671        do l=1,llm
672          do i=1,klon
673            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
674          enddo
675        enddo
676      enddo
677       
678      do l=1,llm
679        do i=1,klon
680          zdufi_omp(i,l)=zdufi(offset+i,l)
681        enddo
682      enddo
683       
684      do l=1,llm
685        do i=1,klon
686          zdvfi_omp(i,l)=zdvfi(offset+i,l)
687        enddo
688      enddo
689       
690      do l=1,llm
691        do i=1,klon
692          zdtfi_omp(i,l)=zdtfi(offset+i,l)
693        enddo
694      enddo
695       
696      do iq=1,nqtot
697        do l=1,llm
698          do i=1,klon
699            zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
700          enddo
701        enddo
702      enddo
703             
704      do i=1,klon
705        zdpsrf_omp(i)=zdpsrf(offset+i)
706      enddo
707
708      do l=1,llm
709        do i=1,klon
710          flxwfi_omp(i,l)=flxwfi(offset+i,l)
711        enddo
712      enddo
713     
714c$OMP BARRIER
715     
716
717!$OMP MASTER
718!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
719!$OMP END MASTER
720      zdt_split=dtphys/nsplit_phys
721      zdufic_omp(:,:)=0.
722      zdvfic_omp(:,:)=0.
723      zdtfic_omp(:,:)=0.
724      zdqfic_omp(:,:,:)=0.
725
726#ifdef CPP_PHYS
727      do isplit=1,nsplit_phys
728
729         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
730         debut_split=debut.and.isplit==1
731         lafin_split=lafin.and.isplit==nsplit_phys
732
733        CALL call_physiq(klon,llm,nqtot,tname,
734     &                   debut_split,lafin_split,
735     &                   jD_cur,jH_cur_split,zdt_split,
736     &                   zplev_omp,zplay_omp,
737     &                   zpk_omp,zphi_omp,zphis_omp,
738     &                   presnivs_omp,
739     &                   zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
740     &                   flxwfi_omp,pducov,
741     &                   zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
742     &                   zdpsrf_omp)
743
744
745         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
746         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
747         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
748         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
749
750         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
751         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
752         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
753         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
754
755      enddo
756
757#endif
758! of #ifdef CPP_PHYS
759
760
761      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
762      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
763      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
764      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
765
766c$OMP BARRIER
767
768      do l=1,llm+1
769        do i=1,klon
770          zplev(offset+i,l)=zplev_omp(i,l)
771        enddo
772      enddo
773         
774       do l=1,llm
775        do i=1,klon 
776          zplay(offset+i,l)=zplay_omp(i,l)
777        enddo
778      enddo
779       
780      do l=1,llm
781        do i=1,klon
782          zphi(offset+i,l)=zphi_omp(i,l)
783        enddo
784      enddo
785       
786
787      do i=1,klon
788        zphis(offset+i)=zphis_omp(i)
789      enddo
790     
791       
792      do l=1,llm
793        presnivs(l)=presnivs_omp(l)
794      enddo
795       
796      do l=1,llm
797        do i=1,klon
798          zufi(offset+i,l)=zufi_omp(i,l)
799        enddo
800      enddo
801       
802      do l=1,llm
803        do i=1,klon
804          zvfi(offset+i,l)=zvfi_omp(i,l)
805        enddo
806      enddo
807       
808      do l=1,llm
809        do i=1,klon
810          ztfi(offset+i,l)=ztfi_omp(i,l)
811        enddo
812      enddo
813       
814      do iq=1,nqtot
815        do l=1,llm
816          do i=1,klon
817            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
818          enddo
819        enddo
820      enddo
821       
822      do l=1,llm
823        do i=1,klon
824          zdufi(offset+i,l)=zdufi_omp(i,l)
825        enddo
826      enddo
827       
828      do l=1,llm
829        do i=1,klon
830          zdvfi(offset+i,l)=zdvfi_omp(i,l)
831        enddo
832      enddo
833       
834      do l=1,llm
835        do i=1,klon
836          zdtfi(offset+i,l)=zdtfi_omp(i,l)
837        enddo
838      enddo
839       
840      do iq=1,nqtot
841        do l=1,llm
842          do i=1,klon
843            zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
844          enddo
845        enddo
846      enddo
847             
848      do i=1,klon
849        zdpsrf(offset+i)=zdpsrf_omp(i)
850      enddo
851     
852
853      klon=klon_mpi
854500   CONTINUE
855c$OMP BARRIER
856
857c$OMP MASTER
858      call stop_timer(timer_physic)
859c$OMP END MASTER
860
861      IF (using_mpi) THEN
862           
863      if (MPI_rank>0) then
864
865c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
866       DO l=1,llm     
867        du_send(1:iim,l)=zdufi(1:iim,l)
868        dv_send(1:iim,l)=zdvfi(1:iim,l)
869       ENDDO
870c$OMP END DO NOWAIT       
871
872c$OMP BARRIER
873#ifdef CPP_MPI
874c$OMP MASTER
875!$OMP CRITICAL (MPI)
876        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
877     &                   COMM_LMDZ,Req(1),ierr)
878        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
879     &                  COMM_LMDZ,Req(2),ierr)
880!$OMP END CRITICAL (MPI)
881c$OMP END MASTER
882#endif
883c$OMP BARRIER
884     
885      endif
886   
887      if (MPI_rank<MPI_Size-1) then
888c$OMP BARRIER
889#ifdef CPP_MPI
890c$OMP MASTER     
891!$OMP CRITICAL (MPI)
892        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
893     &                 COMM_LMDZ,Req(3),ierr)
894        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
895     &                 COMM_LMDZ,Req(4),ierr)
896!$OMP END CRITICAL (MPI)
897c$OMP END MASTER
898#endif
899      endif
900
901c$OMP BARRIER
902
903
904#ifdef CPP_MPI
905c$OMP MASTER   
906!$OMP CRITICAL (MPI)
907      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
908        call MPI_WAITALL(4,Req(1),Status,ierr)
909      else if (MPI_rank>0) then
910        call MPI_WAITALL(2,Req(1),Status,ierr)
911      else if (MPI_rank <MPI_Size-1) then
912        call MPI_WAITALL(2,Req(3),Status,ierr)
913      endif
914!$OMP END CRITICAL (MPI)
915c$OMP END MASTER
916#endif
917
918c$OMP BARRIER     
919
920      ENDIF ! using_mpi
921     
922     
923c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
924      DO l=1,llm
925           
926        zdufi2(1:klon,l)=zdufi(1:klon,l)
927        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
928           
929        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
930        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
931
932        pdhfi(:,jj_begin,l)=0
933        pdqfi(:,jj_begin,l,:)=0
934        pdufi(:,jj_begin,l)=0
935        pdvfi(:,jj_begin,l)=0
936               
937        if (.not. is_south_pole_dyn) then
938          pdhfi(:,jj_end:jj_end+1,l)=0
939          pdqfi(:,jj_end:jj_end+1,l,:)=0
940          pdufi(:,jj_end:jj_end+1,l)=0
941          pdvfi(:,jj_end:jj_end+1,l)=0
942        endif
943     
944       ENDDO
945c$OMP END DO NOWAIT
946
947c$OMP MASTER
948        pdpsfi(:,jj_begin)=0   
949       
950       if (.not. is_south_pole_dyn) then
951         pdpsfi(:,jj_end:jj_end+1)=0
952       endif
953c$OMP END MASTER
954c-----------------------------------------------------------------------
955c   transformation des tendances physiques en tendances dynamiques:
956c   ---------------------------------------------------------------
957
958c  tendance sur la pression :
959c  -----------------------------------
960c      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
961
962c$OMP MASTER
963      kstart=1
964      kend=klon
965
966      if (is_north_pole_dyn) kstart=2
967      if (is_south_pole_dyn)  kend=klon-1
968
969!CDIR ON_ADB(index_i)
970!CDIR ON_ADB(index_j)
971!cdir NODEP
972        do ig0=kstart,kend
973          i=index_i(ig0)
974          j=index_j(ig0)
975          pdpsfi(i,j) = zdpsrf(ig0)
976          if (i==1) pdpsfi(iip1,j) =  zdpsrf(ig0)
977         enddo         
978
979        if (is_north_pole_dyn) then
980            DO i=1,iip1
981              pdpsfi(i,1)    = zdpsrf(1)
982            enddo
983        endif
984       
985        if (is_south_pole_dyn) then
986            DO i=1,iip1
987              pdpsfi(i,jjp1) = zdpsrf(klon)
988            ENDDO
989        endif
990c$OMP END MASTER
991cc$OMP BARRIER
992
993c
994c   62. enthalpie potentielle
995c   ---------------------
996     
997      kstart=1
998      kend=klon
999
1000      if (is_north_pole_dyn) kstart=2
1001      if (is_south_pole_dyn)  kend=klon-1
1002
1003c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1004      DO l=1,llm
1005
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          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1013          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
1014         enddo         
1015
1016        if (is_north_pole_dyn) then
1017            DO i=1,iip1
1018              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
1019            enddo
1020        endif
1021       
1022        if (is_south_pole_dyn) then
1023            DO i=1,iip1
1024              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
1025            ENDDO
1026        endif
1027      ENDDO
1028c$OMP END DO NOWAIT
1029     
1030c   62. humidite specifique
1031c   ---------------------
1032! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1033!      DO iq=1,nqtot
1034!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1035!         DO l=1,llm
1036!!!cdir NODEP
1037!           do ig0=kstart,kend
1038!             i=index_i(ig0)
1039!             j=index_j(ig0)
1040!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1041!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1042!           enddo
1043!           
1044!           if (is_north_pole_dyn) then
1045!             do i=1,iip1
1046!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
1047!             enddo
1048!           endif
1049!           
1050!           if (is_south_pole_dyn) then
1051!             do i=1,iip1
1052!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1053!             enddo
1054!           endif
1055!         ENDDO
1056!c$OMP END DO NOWAIT
1057!      ENDDO
1058
1059c   63. traceurs
1060c   ------------
1061C     initialisation des tendances
1062
1063c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1064      DO l=1,llm
1065        pdqfi(:,jj_begin:jj_end,l,:)=0.
1066      ENDDO
1067c$OMP END DO NOWAIT         
1068
1069C
1070!cdir NODEP
1071      DO iq=1,nqtot
1072         iiq=niadv(iq)
1073c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1074         DO l=1,llm
1075!CDIR ON_ADB(index_i)
1076!CDIR ON_ADB(index_j)
1077!cdir NODEP           
1078             DO ig0=kstart,kend
1079              i=index_i(ig0)
1080              j=index_j(ig0)
1081              pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1082              if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1083            ENDDO
1084           
1085            IF (is_north_pole_dyn) then
1086              DO i=1,iip1
1087                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
1088              ENDDO
1089            ENDIF
1090           
1091            IF (is_south_pole_dyn) then
1092              DO i=1,iip1
1093                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1094              ENDDO
1095            ENDIF
1096           
1097         ENDDO
1098c$OMP END DO NOWAIT         
1099      ENDDO
1100     
1101c   65. champ u:
1102c   ------------
1103c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1104      DO l=1,llm
1105!CDIR ON_ADB(index_i)
1106!CDIR ON_ADB(index_j)
1107!cdir NODEP
1108         do ig0=kstart,kend
1109           i=index_i(ig0)
1110           j=index_j(ig0)
1111           
1112           if (i/=iim) then
1113             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1114           endif
1115           
1116           if (i==1) then
1117              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
1118     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
1119             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1120           endif
1121         
1122         enddo
1123         
1124         if (is_north_pole_dyn) then
1125           DO i=1,iip1
1126            pdufi(i,1,l)    = 0.
1127           ENDDO
1128         endif
1129         
1130         if (is_south_pole_dyn) then
1131           DO i=1,iip1
1132            pdufi(i,jjp1,l) = 0.
1133           ENDDO
1134         endif
1135         
1136      ENDDO
1137c$OMP END DO NOWAIT
1138
1139c   67. champ v:
1140c   ------------
1141
1142      kstart=1
1143      kend=klon
1144
1145      if (is_north_pole_dyn) kstart=2
1146      if (is_south_pole_dyn)  kend=klon-1-iim
1147     
1148c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1149      DO l=1,llm
1150!CDIR ON_ADB(index_i)
1151!CDIR ON_ADB(index_j)
1152!cdir NODEP
1153        do ig0=kstart,kend
1154           i=index_i(ig0)
1155           j=index_j(ig0)
1156           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
1157           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1158     $                                            zdvfi2(ig0+iim,l))
1159     $                                          *cv(i,j)
1160        enddo
1161         
1162      ENDDO
1163c$OMP END DO NOWAIT
1164
1165
1166c   68. champ v pres des poles:
1167c   ---------------------------
1168c      v = U * cos(long) + V * SIN(long)
1169
1170      if (is_north_pole_dyn) then
1171
1172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1173        DO l=1,llm
1174
1175          DO i=1,iim
1176            pdvfi(i,1,l)=
1177     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
1178       
1179            pdvfi(i,1,l)=
1180     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
1181          ENDDO
1182
1183          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
1184
1185        ENDDO
1186c$OMP END DO NOWAIT
1187
1188      endif   
1189     
1190      if (is_south_pole_dyn) then
1191
1192c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1193         DO l=1,llm
1194 
1195           DO i=1,iim
1196              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1197     $        +zdvfi(klon,l)*SIN(rlonv(i))
1198
1199              pdvfi(i,jjm,l)=
1200     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
1201           ENDDO
1202
1203           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1204
1205        ENDDO
1206c$OMP END DO NOWAIT
1207     
1208      endif
1209c-----------------------------------------------------------------------
1210
1211700   CONTINUE
1212 
1213      firstcal = .FALSE.
1214
1215#else
1216      write(lunout,*)
1217     & "calfis_p: for now can only work with parallel physics"
1218      stop
1219#endif
1220! of #ifdef CPP_PHYS
1221#endif
1222! of #ifdef CPP_PARA
1223      END
Note: See TracBrowser for help on using the repository browser.