source: LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F @ 1763

Last change on this file since 1763 was 1749, checked in by Ehouarn Millour, 12 years ago

Added handling of Newtonian (and Shallow Water) modes in dyn3dmem:

  • adapted calfis_loc.F and gcm.F
  • removed unused routines divgrad_p.F and gradiv_p.F
  • adapted iniacademic.F90 and sw_case_williamson.F to work in parallel and renamed them iniacademi_loc.F90 and sw_case_williamson_loc.F
  • fixed bug in exner_milieu_loc.F (filtreg_p form mod_filtreg_p should be used)

EM

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