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