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

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

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