source: LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F @ 4643

Last change on this file since 4643 was 4600, checked in by yann meurdesoif, 17 months ago

Suppress CPP_MPI key usage in source code. MPI wrappers is used to supress missing symbol if the mpi library is not linked

YM

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