source: LMDZ5/trunk/libf/dynphy_lonlat/calfis_loc.F @ 2418

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

Improving the physics/dynamics interface:

  • added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
  • turned "physiq.F90" into module "physiq_mod.F90" for better control of "physiq" arguments. Extracted embeded "gr_fi_ecrit" as self-standing routine (but note that this routine actually only works in serial mode).

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