source: LMDZ5/trunk/libf/dyn3dpar/calfis_p.F @ 2031

Last change on this file since 2031 was 1987, checked in by Ehouarn Millour, 11 years ago

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
EM

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