source: LMDZ6/branches/contrails/libf/dynphy_lonlat/calfis_loc.F90 @ 5428

Last change on this file since 5428 was 5367, checked in by abarral, 3 weeks ago

(WIP) Turn implicit into explicit declarations

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