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

Last change on this file since 1657 was 1657, checked in by Laurent Fairhead, 12 years ago

Phasage de la dynamique parallèle localisée (petite mémoire) avec la branche LMDZ4V5.0-dev (fin de la branche)
Validation effectuée par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modèle donnée paY. Meurdesoif et la version phasée avec la r1399 (fin de la branche LMDZ4V5.0-dev)

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