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