source: LMDZ4/trunk/libf/dyn3dpar/calfis_p.F @ 2206

Last change on this file since 2206 was 1407, checked in by Laurent Fairhead, 14 years ago

Some extraneous prints deleted for optimisation purposes


Des prints supprimés pour optimisation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.0 KB
RevLine 
[630]1!
[1279]2! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fhourdin $
[630]3!
4C
5C
[1146]6      SUBROUTINE calfis_p(lafin,
[1279]7     $                  jD_cur, jH_cur,
[630]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)
[1279]29#ifdef CPP_EARTH
30! Ehouarn: For now, calfis_p needs Earth physics
[630]31c
32c    Auteur :  P. Le Van, F. Hourdin
33c   .........
34      USE dimphy
[774]35      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
[1403]36      USE mod_interface_dyn_phys
37      USE IOPHY
38#endif
[1000]39      USE parallel, ONLY : omp_chunk, using_mpi
[630]40      USE Write_Field
41      Use Write_field_p
42      USE Times
[1146]43      USE infotrac
[1403]44      USE control_mod
[1146]45
[630]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
[1146]105      INTEGER ngridmx
[630]106      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
107
108#include "comconst.h"
109#include "comvert.h"
110#include "comgeom2.h"
[1403]111#include "iniprint.h"
[1000]112#ifdef CPP_MPI
[630]113      include 'mpif.h'
[1000]114#endif
[630]115c    Arguments :
116c    -----------
117      LOGICAL  lafin
[1403]118!      REAL heure
119      REAL, intent(in):: jD_cur, jH_cur
[630]120      REAL pvcov(iip1,jjm,llm)
121      REAL pucov(iip1,jjp1,llm)
122      REAL pteta(iip1,jjp1,llm)
123      REAL pmasse(iip1,jjp1,llm)
[1146]124      REAL pq(iip1,jjp1,llm,nqtot)
[630]125      REAL pphis(iip1,jjp1)
126      REAL pphi(iip1,jjp1,llm)
127c
128      REAL pdvcov(iip1,jjm,llm)
129      REAL pducov(iip1,jjp1,llm)
130      REAL pdteta(iip1,jjp1,llm)
[1146]131      REAL pdq(iip1,jjp1,llm,nqtot)
[1403]132      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
[630]133c
134      REAL pps(iip1,jjp1)
135      REAL pp(iip1,jjp1,llmp1)
136      REAL ppk(iip1,jjp1,llm)
137c
138      REAL pdvfi(iip1,jjm,llm)
139      REAL pdufi(iip1,jjp1,llm)
140      REAL pdhfi(iip1,jjp1,llm)
[1146]141      REAL pdqfi(iip1,jjp1,llm,nqtot)
[630]142      REAL pdpsfi(iip1,jjp1)
143
144      INTEGER        longcles
145      PARAMETER    ( longcles = 20 )
146      REAL clesphy0( longcles )
147
[1403]148#ifdef CPP_EARTH
[630]149c    Local variables :
150c    -----------------
151
152      INTEGER i,j,l,ig0,ig,iq,iiq
[764]153      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
154      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
155      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
[630]156c
[764]157      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
158      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
[630]159c
[764]160      REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
161      REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
[630]162c
[764]163      REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
164      REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
165      REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
[985]166      REAL,SAVE,ALLOCATABLE ::  flxwfi(:,:)     ! Flux de masse verticale sur la grille physiq
167
[630]168c
[764]169      REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
170      REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
171      REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
172      REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
173      REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
174      REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
175      REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
176      REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
177      REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
178      REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
179      REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
180      REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
181      REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
182      REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
[985]183      REAL,SAVE,ALLOCATABLE ::  flxwfi_omp(:,:)     ! Flux de masse verticale sur la grille physiq
[764]184
[1403]185!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186! Introduction du splitting (FH)
187! Question pour Yann :
188! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
189! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
190! soit allocatable (plutot par exemple que de passer une dimension
191! dépendant du process en argument des routines) et que, du coup,
192! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
193! Tu confirmes ?
194! J'ai suivi le même principe pour les zdufic_omp
195! Mais c'est surement bien que tu controles.
196!
197
198      REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
199      REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
200      REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
201      REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
202      REAL jH_cur_split,zdt_split
203      LOGICAL debut_split,lafin_split
204      INTEGER isplit
205!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
206
[764]207c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
208c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
[985]209c$OMP+                 zqfi_omp,zdufi_omp,zdvfi_omp,
[1403]210c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
211c$OMP+                 zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)       
[764]212
213      LOGICAL,SAVE :: first_omp=.true.
214c$OMP THREADPRIVATE(first_omp)
215     
[630]216      REAL zsin(iim),zcos(iim),z1(iim)
217      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
218      REAL unskap, pksurcp
[764]219c
220cIM diagnostique PVteta, Amip2
221      INTEGER ntetaSTD
222      PARAMETER(ntetaSTD=3)
223      REAL rtetaSTD(ntetaSTD)
224      DATA rtetaSTD/350., 380., 405./
225      REAL PVteta(klon,ntetaSTD)
226     
[630]227     
228      REAL SSUM
229
230      LOGICAL firstcal, debut
231      DATA firstcal/.true./
232      SAVE firstcal,debut
[764]233c$OMP THREADPRIVATE(firstcal,debut)
[630]234     
[764]235      REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
[630]236      INTEGER :: ierr
[1000]237#ifdef CPP_MPI
[630]238      INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
[1000]239#else
240      INTEGER,dimension(1,4) :: Status
241#endif
[630]242      INTEGER, dimension(4) :: Req
[764]243      REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
244      integer :: k,kstart,kend
245      INTEGER :: offset 
[630]246c
247c-----------------------------------------------------------------------
248c
249c    1. Initialisations :
250c    --------------------
251c
252
[764]253      klon=klon_mpi
254     
255      PVteta(:,:)=0.
256           
[1279]257c
258      IF ( firstcal )  THEN
259        debut = .TRUE.
260        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
[1403]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
[630]266         STOP
[1279]267        ENDIF
[764]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))
[1146]273      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
[764]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))
[1146]277      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
[764]278      ALLOCATE(zdpsrf(klon))
279      ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
[985]280      ALLOCATE(flxwfi(klon,llm))
[764]281c$OMP END MASTER
282c$OMP BARRIER     
[630]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
[764]296c$OMP MASTER
[630]297      call start_timer(timer_physic)
[764]298c$OMP END MASTER
299
300c$OMP MASTER             
[1279]301!CDIR ON_ADB(index_i)
302!CDIR ON_ADB(index_j)
[630]303      do ig0=1,klon
[774]304        i=index_i(ig0)
305        j=index_j(ig0)
[630]306        zpsrf(ig0)=pps(i,j)
307      enddo
[764]308c$OMP END MASTER
[630]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
[961]322c      print *,omp_rank,'klon--->',klon
[764]323c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]324      DO l = 1, llmp1
[1279]325!CDIR ON_ADB(index_i)
326!CDIR ON_ADB(index_j)
[630]327        do ig0=1,klon
[774]328          i=index_i(ig0)
329          j=index_j(ig0)
[630]330          zplev( ig0,l ) = pp(i,j,l)
331        enddo
332      ENDDO
[764]333c$OMP END DO NOWAIT
[630]334c
335c
336
337c   43. temperature naturelle (en K) et pressions milieux couches .
338c   ---------------------------------------------------------------
[764]339c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]340      DO l=1,llm
[1279]341!CDIR ON_ADB(index_i)
342!CDIR ON_ADB(index_j)
[630]343        do ig0=1,klon
[774]344          i=index_i(ig0)
345          j=index_j(ig0)
[630]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
[764]352c$OMP END DO NOWAIT
[630]353
354c   43.bis traceurs
355c   ---------------
356c
357
[1146]358      DO iq=1,nqtot
[630]359         iiq=niadv(iq)
[764]360c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]361         DO l=1,llm
[1279]362!CDIR ON_ADB(index_i)
363!CDIR ON_ADB(index_j)
[630]364           do ig0=1,klon
[774]365             i=index_i(ig0)
366             j=index_j(ig0)
[630]367             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
368           enddo
369         ENDDO
[764]370c$OMP END DO NOWAIT     
[630]371      ENDDO
372
373
374c   Geopotentiel calcule par rapport a la surface locale:
375c   -----------------------------------------------------
376
377      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
[764]378
[630]379      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
[1146]380
[764]381c$OMP BARRIER
[630]382
[764]383c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]384      DO l=1,llm
385         DO ig=1,klon
386           zphi(ig,l)=zphi(ig,l)-zphis(ig)
387         ENDDO
388      ENDDO
[960]389c$OMP END DO NOWAIT
390     
[630]391
392c
393c   45. champ u:
394c   ------------
395
396      kstart=1
397      kend=klon
398     
[774]399      if (is_north_pole) kstart=2
400      if (is_south_pole) kend=klon-1
[630]401     
[764]402c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]403      DO l=1,llm
[1279]404!CDIR ON_ADB(index_i)
405!CDIR ON_ADB(index_j)
406!CDIR SPARSE
[630]407        do ig0=kstart,kend
[774]408          i=index_i(ig0)
409          j=index_j(ig0)
[630]410          if (i==1) then
411            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
412     $                         + pucov(1,j,l)/cu(1,j) )
413          else
414            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j)
415     $                       + pucov(i,j,l)/cu(i,j) )
416          endif
417        enddo
418      ENDDO
[764]419c$OMP END DO NOWAIT
[1279]420
[630]421c   46.champ v:
422c   -----------
[1279]423
[764]424c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]425      DO l=1,llm
[1279]426!CDIR ON_ADB(index_i)
427!CDIR ON_ADB(index_j)
[630]428        DO ig0=kstart,kend
[774]429          i=index_i(ig0)
430          j=index_j(ig0)
[630]431          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1)
432     $                       + pvcov(i,j,l)/cv(i,j) )
433   
434         ENDDO
435      ENDDO
[764]436c$OMP END DO NOWAIT
[630]437
438c   47. champs de vents aux pole nord   
439c   ------------------------------
440c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
441c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
442
[774]443      if (is_north_pole) then
[764]444c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]445        DO l=1,llm
446
447           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
448           DO i=2,iim
449              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
450           ENDDO
451 
452           DO i=1,iim
453              zcos(i)   = COS(rlonv(i))*z1(i)
454              zsin(i)   = SIN(rlonv(i))*z1(i)
455           ENDDO
456 
457           zufi(1,l)  = SSUM(iim,zcos,1)/pi
458           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
459 
460        ENDDO
[764]461c$OMP END DO NOWAIT     
[630]462      endif
463
464
465c   48. champs de vents aux pole sud:
466c   ---------------------------------
467c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
468c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
469
[774]470      if (is_south_pole) then
[764]471c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]472        DO l=1,llm
473 
474         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
475           DO i=2,iim
[1279]476             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
[630]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(klon,l)  = SSUM(iim,zcos,1)/pi
485           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
486        ENDDO
[764]487c$OMP END DO NOWAIT       
[630]488      endif
489
490
[774]491      IF (is_sequential) THEN
[764]492c
493cIM calcul PV a teta=350, 380, 405K
494        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
495     $           ztfi,zplay,zplev,
496     $           ntetaSTD,rtetaSTD,PVteta)
497c
498      ENDIF
[960]499
500c On change de grille, dynamique vers physiq, pour le flux de masse verticale
[630]501      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
502
503c-----------------------------------------------------------------------
504c   Appel de la physique:
505c   ---------------------
506
507
[764]508c$OMP BARRIER
509      if (first_omp) then
[774]510        klon=klon_omp
[764]511
512        allocate(zplev_omp(klon,llm+1))
513        allocate(zplay_omp(klon,llm))
514        allocate(zphi_omp(klon,llm))
515        allocate(zphis_omp(klon))
516        allocate(presnivs_omp(llm))
517        allocate(zufi_omp(klon,llm))
518        allocate(zvfi_omp(klon,llm))
519        allocate(ztfi_omp(klon,llm))
[1146]520        allocate(zqfi_omp(klon,llm,nqtot))
[764]521        allocate(zdufi_omp(klon,llm))
522        allocate(zdvfi_omp(klon,llm))
523        allocate(zdtfi_omp(klon,llm))
[1146]524        allocate(zdqfi_omp(klon,llm,nqtot))
[1403]525        allocate(zdufic_omp(klon,llm))
526        allocate(zdvfic_omp(klon,llm))
527        allocate(zdtfic_omp(klon,llm))
528        allocate(zdqfic_omp(klon,llm,nqtot))
[764]529        allocate(zdpsrf_omp(klon))
[985]530        allocate(flxwfi_omp(klon,llm))
[764]531        first_omp=.false.
532      endif
533       
534           
[774]535      klon=klon_omp
536      offset=klon_omp_begin-1
[764]537     
538      do l=1,llm+1
539        do i=1,klon
540          zplev_omp(i,l)=zplev(offset+i,l)
541        enddo
542      enddo
543         
544       do l=1,llm
545        do i=1,klon 
546          zplay_omp(i,l)=zplay(offset+i,l)
547        enddo
548      enddo
549       
550      do l=1,llm
551        do i=1,klon
552          zphi_omp(i,l)=zphi(offset+i,l)
553        enddo
554      enddo
555       
556      do i=1,klon
557        zphis_omp(i)=zphis(offset+i)
558      enddo
559     
560       
561      do l=1,llm
562        presnivs_omp(l)=presnivs(l)
563      enddo
564       
565      do l=1,llm
566        do i=1,klon
567          zufi_omp(i,l)=zufi(offset+i,l)
568        enddo
569      enddo
570       
571      do l=1,llm
572        do i=1,klon
573          zvfi_omp(i,l)=zvfi(offset+i,l)
574        enddo
575      enddo
576       
577      do l=1,llm
578        do i=1,klon
579          ztfi_omp(i,l)=ztfi(offset+i,l)
580        enddo
581      enddo
582       
[1146]583      do iq=1,nqtot
[764]584        do l=1,llm
585          do i=1,klon
586            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
587          enddo
588        enddo
589      enddo
590       
591      do l=1,llm
592        do i=1,klon
593          zdufi_omp(i,l)=zdufi(offset+i,l)
594        enddo
595      enddo
596       
597      do l=1,llm
598        do i=1,klon
599          zdvfi_omp(i,l)=zdvfi(offset+i,l)
600        enddo
601      enddo
602       
603      do l=1,llm
604        do i=1,klon
605          zdtfi_omp(i,l)=zdtfi(offset+i,l)
606        enddo
607      enddo
608       
[1146]609      do iq=1,nqtot
[764]610        do l=1,llm
611          do i=1,klon
612            zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
613          enddo
614        enddo
615      enddo
616       
617      do i=1,klon
618        zdpsrf_omp(i)=zdpsrf(offset+i)
619      enddo
[985]620
621      do l=1,llm
622        do i=1,klon
623          flxwfi_omp(i,l)=flxwfi(offset+i,l)
624        enddo
625      enddo
[764]626     
627c$OMP BARRIER
628     
[1279]629      if (planet_type=="earth") then
630#ifdef CPP_EARTH
[1403]631
632!$OMP MASTER
[1407]633!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
[1403]634!$OMP END MASTER
635      zdt_split=dtphys/nsplit_phys
636      zdufic_omp(:,:)=0.
637      zdvfic_omp(:,:)=0.
638      zdtfic_omp(:,:)=0.
639      zdqfic_omp(:,:,:)=0.
640
641      do isplit=1,nsplit_phys
642
643         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
644         debut_split=debut.and.isplit==1
645         lafin_split=lafin.and.isplit==nsplit_phys
646
647
[630]648      CALL physiq (klon,
649     .             llm,
[1403]650     .             debut_split,
651     .             lafin_split,
[1279]652     .             jD_cur,
[1403]653     .             jH_cur_split,
654     .             zdt_split,
[764]655     .             zplev_omp,
656     .             zplay_omp,
657     .             zphi_omp,
658     .             zphis_omp,
659     .             presnivs_omp,
[630]660     .             clesphy0,
[764]661     .             zufi_omp,
662     .             zvfi_omp,
663     .             ztfi_omp,
664     .             zqfi_omp,
[960]665c#ifdef INCA
[985]666     .             flxwfi_omp,
[960]667c#endif
[764]668     .             zdufi_omp,
669     .             zdvfi_omp,
670     .             zdtfi_omp,
671     .             zdqfi_omp,
672     .             zdpsrf_omp,
673cIM diagnostique PVteta, Amip2         
674     .             pducov,
675     .             PVteta)
[1403]676
677         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
678         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
679         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
680         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
681
682         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
683         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
684         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
685         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
686
687      enddo
688
689      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
690      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
691      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
692      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
693
[1279]694#endif
695      endif !of if (planet_type=="earth")
[764]696c$OMP BARRIER
697
698      do l=1,llm+1
699        do i=1,klon
700          zplev(offset+i,l)=zplev_omp(i,l)
701        enddo
702      enddo
703         
704       do l=1,llm
705        do i=1,klon 
706          zplay(offset+i,l)=zplay_omp(i,l)
707        enddo
708      enddo
709       
710      do l=1,llm
711        do i=1,klon
712          zphi(offset+i,l)=zphi_omp(i,l)
713        enddo
714      enddo
715       
716
717      do i=1,klon
718        zphis(offset+i)=zphis_omp(i)
719      enddo
720     
721       
722      do l=1,llm
723        presnivs(l)=presnivs_omp(l)
724      enddo
725       
726      do l=1,llm
727        do i=1,klon
728          zufi(offset+i,l)=zufi_omp(i,l)
729        enddo
730      enddo
731       
732      do l=1,llm
733        do i=1,klon
734          zvfi(offset+i,l)=zvfi_omp(i,l)
735        enddo
736      enddo
737       
738      do l=1,llm
739        do i=1,klon
740          ztfi(offset+i,l)=ztfi_omp(i,l)
741        enddo
742      enddo
743       
[1146]744      do iq=1,nqtot
[764]745        do l=1,llm
746          do i=1,klon
747            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
748          enddo
749        enddo
750      enddo
751       
752      do l=1,llm
753        do i=1,klon
754          zdufi(offset+i,l)=zdufi_omp(i,l)
755        enddo
756      enddo
757       
758      do l=1,llm
759        do i=1,klon
760          zdvfi(offset+i,l)=zdvfi_omp(i,l)
761        enddo
762      enddo
763       
764      do l=1,llm
765        do i=1,klon
766          zdtfi(offset+i,l)=zdtfi_omp(i,l)
767        enddo
768      enddo
769       
[1146]770      do iq=1,nqtot
[764]771        do l=1,llm
772          do i=1,klon
773            zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
774          enddo
775        enddo
776      enddo
777       
778      do i=1,klon
779        zdpsrf(offset+i)=zdpsrf_omp(i)
780      enddo
781     
782
783      klon=klon_mpi
[630]784500   CONTINUE
[764]785c$OMP BARRIER
[630]786
[764]787c$OMP MASTER
[630]788      call stop_timer(timer_physic)
[764]789c$OMP END MASTER
[1000]790
791      IF (using_mpi) THEN
792           
[630]793      if (MPI_rank>0) then
[764]794
795c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
796       DO l=1,llm     
797        du_send(1:iim,l)=zdufi(1:iim,l)
798        dv_send(1:iim,l)=zdvfi(1:iim,l)
799       ENDDO
800c$OMP END DO NOWAIT       
801
802c$OMP BARRIER
[1000]803#ifdef CPP_MPI
[764]804c$OMP MASTER
[985]805!$OMP CRITICAL (MPI)
[630]806        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
[764]807     &                   COMM_LMDZ,Req(1),ierr)
[630]808        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
[764]809     &                  COMM_LMDZ,Req(2),ierr)
[985]810!$OMP END CRITICAL (MPI)
[764]811c$OMP END MASTER
[1000]812#endif
[764]813c$OMP BARRIER
[630]814     
815      endif
816   
817      if (MPI_rank<MPI_Size-1) then
[764]818c$OMP BARRIER
[1000]819#ifdef CPP_MPI
[764]820c$OMP MASTER     
[985]821!$OMP CRITICAL (MPI)
[630]822        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
[764]823     &                 COMM_LMDZ,Req(3),ierr)
[630]824        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
[764]825     &                 COMM_LMDZ,Req(4),ierr)
[985]826!$OMP END CRITICAL (MPI)
[764]827c$OMP END MASTER
[1000]828#endif
[630]829      endif
[764]830
831c$OMP BARRIER
[1000]832
833
834#ifdef CPP_MPI
[764]835c$OMP MASTER   
[985]836!$OMP CRITICAL (MPI)
[630]837      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
838        call MPI_WAITALL(4,Req(1),Status,ierr)
839      else if (MPI_rank>0) then
840        call MPI_WAITALL(2,Req(1),Status,ierr)
841      else if (MPI_rank <MPI_Size-1) then
842        call MPI_WAITALL(2,Req(3),Status,ierr)
843      endif
[985]844!$OMP END CRITICAL (MPI)
[764]845c$OMP END MASTER
[1000]846#endif
847
[764]848c$OMP BARRIER     
849
[1000]850      ENDIF ! using_mpi
851     
852     
[764]853c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
854      DO l=1,llm
855           
856        zdufi2(1:klon,l)=zdufi(1:klon,l)
857        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
858           
859        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
860        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
861 
[774]862         pdhfi(:,jj_begin,l)=0
863         pdqfi(:,jj_begin,l,:)=0
864         pdufi(:,jj_begin,l)=0
865         pdvfi(:,jj_begin,l)=0
[764]866         
[774]867         if (.not. is_south_pole) then
868           pdhfi(:,jj_end,l)=0
869           pdqfi(:,jj_end,l,:)=0
870           pdufi(:,jj_end,l)=0
871           pdvfi(:,jj_end,l)=0
[764]872         endif
[630]873     
[764]874       ENDDO
875c$OMP END DO NOWAIT
[630]876
[764]877c$OMP MASTER
[774]878       pdpsfi(:,jj_begin)=0   
879       if (.not. is_south_pole) then
880         pdpsfi(:,jj_end)=0
[630]881       endif
[764]882c$OMP END MASTER
[630]883c-----------------------------------------------------------------------
884c   transformation des tendances physiques en tendances dynamiques:
885c   ---------------------------------------------------------------
886
887c  tendance sur la pression :
888c  -----------------------------------
889      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
890c
891c   62. enthalpie potentielle
892c   ---------------------
893     
894      kstart=1
895      kend=klon
896
[774]897      if (is_north_pole) kstart=2
898      if (is_south_pole)  kend=klon-1
[630]899
[764]900c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]901      DO l=1,llm
902
[1279]903!CDIR ON_ADB(index_i)
904!CDIR ON_ADB(index_j)
905!cdir NODEP
[630]906        do ig0=kstart,kend
[774]907          i=index_i(ig0)
908          j=index_j(ig0)
[630]909          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
910          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
911         enddo         
912
[774]913        if (is_north_pole) then
[630]914            DO i=1,iip1
915              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
916            enddo
917        endif
918       
[774]919        if (is_south_pole) then
[630]920            DO i=1,iip1
921              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
922            ENDDO
923        endif
924      ENDDO
[764]925c$OMP END DO NOWAIT
[630]926     
927c   62. humidite specifique
928c   ---------------------
[1279]929! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
930!      DO iq=1,nqtot
931!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
932!         DO l=1,llm
933!!!cdir NODEP
934!           do ig0=kstart,kend
935!             i=index_i(ig0)
936!             j=index_j(ig0)
937!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
938!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
939!           enddo
940!           
941!           if (is_north_pole) then
942!             do i=1,iip1
943!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
944!             enddo
945!           endif
946!           
947!           if (is_south_pole) then
948!             do i=1,iip1
949!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
950!             enddo
951!           endif
952!         ENDDO
953!c$OMP END DO NOWAIT
954!      ENDDO
[630]955
956c   63. traceurs
957c   ------------
958C     initialisation des tendances
[764]959
960c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
961      DO l=1,llm
962        pdqfi(:,:,l,:)=0.
963      ENDDO
964c$OMP END DO NOWAIT     
965
[630]966C
[1279]967!cdir NODEP
[1146]968      DO iq=1,nqtot
[630]969         iiq=niadv(iq)
[764]970c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]971         DO l=1,llm
[1279]972!CDIR ON_ADB(index_i)
973!CDIR ON_ADB(index_j)
974!cdir NODEP           
[630]975             DO ig0=kstart,kend
[774]976              i=index_i(ig0)
977              j=index_j(ig0)
[630]978              pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
979              if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
980            ENDDO
981           
[774]982            IF (is_north_pole) then
[630]983              DO i=1,iip1
984                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
985              ENDDO
986            ENDIF
987           
[774]988            IF (is_south_pole) then
[630]989              DO i=1,iip1
990                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
991              ENDDO
992            ENDIF
993           
994         ENDDO
[764]995c$OMP END DO NOWAIT     
[630]996      ENDDO
997     
998c   65. champ u:
999c   ------------
[764]1000c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]1001      DO l=1,llm
[1279]1002!CDIR ON_ADB(index_i)
1003!CDIR ON_ADB(index_j)
1004!cdir NODEP
[630]1005         do ig0=kstart,kend
[774]1006           i=index_i(ig0)
1007           j=index_j(ig0)
[630]1008           
1009           if (i/=iim) then
1010             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
1011           endif
1012           
1013           if (i==1) then
1014              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
1015     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
[764]1016             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
[630]1017           endif
1018         
1019         enddo
1020         
[774]1021         if (is_north_pole) then
[630]1022           DO i=1,iip1
1023            pdufi(i,1,l)    = 0.
1024           ENDDO
1025         endif
1026         
[774]1027         if (is_south_pole) then
[630]1028           DO i=1,iip1
1029            pdufi(i,jjp1,l) = 0.
1030           ENDDO
1031         endif
1032         
1033      ENDDO
[764]1034c$OMP END DO NOWAIT
[630]1035
1036c   67. champ v:
1037c   ------------
1038
1039      kstart=1
1040      kend=klon
1041
[774]1042      if (is_north_pole) kstart=2
1043      if (is_south_pole)  kend=klon-1-iim
[630]1044     
[764]1045c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]1046      DO l=1,llm
[1279]1047!CDIR ON_ADB(index_i)
1048!CDIR ON_ADB(index_j)
1049!cdir NODEP
[630]1050        do ig0=kstart,kend
[774]1051           i=index_i(ig0)
1052           j=index_j(ig0)
[630]1053           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
1054           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1055     $                                      zdvfi2(ig0+iim,l))
1056     $                                    *cv(i,j)
1057        enddo
1058         
1059      ENDDO
[764]1060c$OMP END DO NOWAIT
[630]1061
1062
1063c   68. champ v pres des poles:
1064c   ---------------------------
1065c      v = U * cos(long) + V * SIN(long)
1066
[774]1067      if (is_north_pole) then
[764]1068
1069c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]1070        DO l=1,llm
1071
1072          DO i=1,iim
1073            pdvfi(i,1,l)=
1074     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
1075       
1076            pdvfi(i,1,l)=
1077     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
1078          ENDDO
1079
1080          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
1081
1082        ENDDO
[764]1083c$OMP END DO NOWAIT
[630]1084
1085      endif   
1086     
[774]1087      if (is_south_pole) then
[764]1088
1089c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1090         DO l=1,llm
[630]1091 
1092           DO i=1,iim
1093              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1094     $        +zdvfi(klon,l)*SIN(rlonv(i))
1095
1096              pdvfi(i,jjm,l)=
1097     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
1098           ENDDO
1099
1100           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1101
1102        ENDDO
[764]1103c$OMP END DO NOWAIT
[630]1104     
1105      endif
1106c-----------------------------------------------------------------------
1107
1108700   CONTINUE
1109 
1110      firstcal = .FALSE.
1111
[1279]1112#else
[1403]1113      write(lunout,*)
1114     & "calfis_p: for now can only work with parallel physics"
[1279]1115      stop
1116#endif
1117! of #ifdef CPP_EARTH
[630]1118      RETURN
1119      END
Note: See TracBrowser for help on using the repository browser.