Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/calfis_p.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 49.2 KB

Line
1	!
2	! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fairhead $
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	$ pdufi,
24	$ pdvfi,
25	$ pdhfi,
26	$ pdqfi,
27	$ pdpsfi)
28
29	! Ehouarn: if using (parallelized) physics
30	USE dimphy
31	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
32	USE mod_interface_dyn_phys
33	! USE IOPHY
34
35	USE parallel_lmdz, ONLY : omp_chunk, using_mpi, AllGather_Field
36	USE Write_Field
37	Use Write_field_p
38	USE Times
39	USE infotrac, ONLY: nqtot, niadv, tname
40	USE control_mod, ONLY: planet_type, nsplit_phys
41	USE cpdet_mod, only: tpot2t_p, t2tpot_p
42
43	! used only for zonal averages
44	USE moyzon_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 (K/s)
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	!-----------------------------------------------------------------------
102	! INCLUDE 'dimensions.h'
103	!
104	! dimensions.h contient les dimensions du modele
105	! ndm est tel que iim=2**ndm
106	!-----------------------------------------------------------------------
107
108	INTEGER iim,jjm,llm,ndm
109
110	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
111
112	!-----------------------------------------------------------------------
113	!
114	! $Header$
115	!
116	!
117	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
118	! veillez n'utiliser que des ! pour les commentaires
119	! et bien positionner les & des lignes de continuation
120	! (les placer en colonne 6 et en colonne 73)
121	!
122	!
123	!-----------------------------------------------------------------------
124	! INCLUDE 'paramet.h'
125
126	INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
127	INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
128	INTEGER ijmllm,mvar
129	INTEGER jcfil,jcfllm
130
131	PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
132	& ,jjp1=jjm+1-1/jjm)
133	PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
134	PARAMETER( kftd = iim/2 -ndm )
135	PARAMETER( ip1jm = iip1jjm, ip1jmp1= iip1jjp1 )
136	PARAMETER( ip1jmi1= ip1jm - iip1 )
137	PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
138	PARAMETER( mvar= ip1jmp1( 2llm+1) + ijmllm )
139	PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
140
141	!-----------------------------------------------------------------------
142	!
143	! $Id: temps.h 1577 2011-10-20 15:06:47Z fairhead $
144	!
145	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
146	! veillez n'utiliser que des ! pour les commentaires
147	! et bien positionner les & des lignes de continuation
148	! (les placer en colonne 6 et en colonne 73)
149	!
150	!
151	! jD_ref = jour julien de la date de reference (lancement de l'experience)
152	! hD_ref = "heure" julienne de la date de reference
153	!-----------------------------------------------------------------------
154	! INCLUDE 'temps.h'
155
156	COMMON/temps_r/dt,jD_ref,jH_ref,start_time,hour_ini
157	COMMON/temps_i/day_ini,day_end,annee_ref,day_ref, &
158	& itau_dyn,itau_phy,itaufin
159	COMMON/temps_c/calend
160
161
162	INTEGER itaufin ! total number of dynamical steps for the run
163	INTEGER itau_dyn, itau_phy
164	INTEGER day_ini ! initial day # of simulation sequence
165	INTEGER day_end ! final day # ; i.e. day # when this simulation ends
166	INTEGER annee_ref
167	INTEGER day_ref
168	REAL dt ! (dynamics) time step (changes if doing Matsuno or LF step)
169	REAL jD_ref, jH_ref, start_time
170	CHARACTER (len=10) :: calend
171
172	! Additionnal Mars stuff:
173	real hour_ini ! initial fraction of day of simulation sequence (0=<hour_ini<1)
174
175	!-----------------------------------------------------------------------
176	!
177	! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
178	!
179	!
180	! NB: keep items of different kinds in seperate common blocs to avoid
181	! "misaligned commons" issues
182	!-----------------------------------------------------------------------
183	! INCLUDE 'logic.h'
184
185	COMMON/logicl/ purmats,forward,leapf,apphys, &
186	& statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus &
187	& ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile &
188	& ,ok_limit,ok_etat0,hybrid &
189	& ,moyzon_mu,moyzon_ch
190
191	COMMON/logici/ iflag_phys,iflag_trac
192
193	LOGICAL purmats,forward,leapf,apphys,statcl,conser, &
194	& apdiss,apdelq,saison,ecripar,fxyhypb,ysinus &
195	& ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile &
196	& ,ok_limit,ok_etat0
197	logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
198	! (only used if disvert_type==2)
199	logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
200
201	integer iflag_phys,iflag_trac
202	!$OMP THREADPRIVATE(/logicl/)
203	!$OMP THREADPRIVATE(/logici/)
204	!-----------------------------------------------------------------------
205
206	INTEGER ngridmx
207	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
208
209	!
210	! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
211	!
212	!-----------------------------------------------------------------------
213	! INCLUDE comconst.h
214
215	COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, &
216	& iflag_top_bound,mode_top_bound
217	COMMON/comconstr/dtvr,daysec, &
218	& pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg &
219	& ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta &
220	& ,dissip_pupstart ,tau_top_bound, &
221	& daylen,molmass, ihf
222	COMMON/cpdetvenus/cpp,nu_venus,t0_venus
223
224	INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
225	REAL dtvr ! dynamical time step (in s)
226	REAL daysec !length (in s) of a standard day
227	REAL pi ! something like 3.14159....
228	REAL dtphys ! (s) time step for the physics
229	REAL dtdiss ! (s) time step for the dissipation
230	REAL rad ! (m) radius of the planet
231	REAL r ! Reduced Gas constant r=R/mu
232	! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
233	REAL cpp ! Cp
234	REAL kappa ! kappa=R/Cp
235	REAL cotot
236	REAL unsim ! = 1./iim
237	REAL g ! (m/s2) gravity
238	REAL omeg ! (rad/s) rotation rate of the planet
239	! Dissipation factors, for Earth model:
240	REAL dissip_factz,dissip_zref !dissip_deltaz
241	! Dissipation factors, for other planets:
242	REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
243	REAL dissip_pupstart
244	INTEGER iflag_top_bound,mode_top_bound
245	REAL tau_top_bound
246	REAL daylen ! length of solar day, in 'standard' day length
247	REAL molmass ! (g/mol) molar mass of the atmosphere
248
249	REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
250	REAL ihf ! (W/m2) intrinsic heat flux for giant planets
251
252
253	!-----------------------------------------------------------------------
254	!
255	! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
256	!
257	!-----------------------------------------------------------------------
258	! INCLUDE 'comvert.h'
259
260	COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), &
261	& pa,preff,nivsigs(llm),nivsig(llm+1), &
262	& aps(llm),bps(llm),scaleheight,pseudoalt(llm)
263
264	common/comverti/disvert_type, pressure_exner
265
266	real ap ! hybrid pressure contribution at interlayers
267	real bp ! hybrid sigma contribution at interlayer
268	real presnivs ! (reference) pressure at mid-layers
269	real dpres
270	real pa ! reference pressure (Pa) at which hybrid coordinates
271	! become purely pressure
272	real preff ! reference surface pressure (Pa)
273	real nivsigs
274	real nivsig
275	real aps ! hybrid pressure contribution at mid-layers
276	real bps ! hybrid sigma contribution at mid-layers
277	real scaleheight ! atmospheric (reference) scale height (km)
278	real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
279	! preff and scaleheight
280
281	integer disvert_type ! type of vertical discretization:
282	! 1: Earth (default for planet_type==earth),
283	! automatic generation
284	! 2: Planets (default for planet_type!=earth),
285	! using 'z2sig.def' (or 'esasig.def) file
286
287	logical pressure_exner
288	! compute pressure inside layers using Exner function, else use mean
289	! of pressure values at interfaces
290
291	!-----------------------------------------------------------------------
292	!
293	! $Header$
294	!
295	!CDK comgeom2
296	COMMON/comgeom/ &
297	& cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm) , &
298	& aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1) , &
299	& airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols , &
300	& unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm) , &
301	& aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1) , &
302	& aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1) , &
303	& alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1) , &
304	& alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1) , &
305	& fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm), &
306	& rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm) , &
307	& cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1) , &
308	& cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
309	& cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
310	& unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2 , &
311	& unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1) , &
312	& unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm) &
313	& , xprimu(iip1),xprimv(iip1)
314
315
316	REAL &
317	& cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
318	& ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4 , &
319	& alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
320	& fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2 , &
321	& cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1 , &
322	& unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2 , &
323	& unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu , &
324	& cusurcvu,xprimu,xprimv
325	!
326	! $Header$
327	!
328	!
329	! gestion des impressions de sorties et de débogage
330	! lunout: unité du fichier dans lequel se font les sorties
331	! (par defaut 6, la sortie standard)
332	! prt_level: niveau d'impression souhaité (0 = minimum)
333	!
334	INTEGER lunout, prt_level
335	COMMON /comprint/ lunout, prt_level
336	include 'mpif.h'
337	c Arguments :
338	c -----------
339	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
340	REAL,INTENT(IN) :: jD_cur, jH_cur
341	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
342	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
343	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
344	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
345	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
346	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
347	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
348
349	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
350	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
351	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
352	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
353	! qui lui meme ne sert a rien dans la routine telle qu'elle est
354	! ecrite, et que j'ai donc commente....
355	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
356	! NB: pdq is only used to compute pcvgq which is in fact not used...
357
358	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
359	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
360	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
361	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on dynamics grid
362
363	! tendencies (in */s) from the physics
364	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
365	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
366	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
367	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
368	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
369
370	c Local variables :
371	c -----------------
372
373	INTEGER i,j,l,ig0,ig,iq,iiq
374	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
375	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
376	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
377
378	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
379	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
380	! ADAPTATION GCM POUR CP(T)
381	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
382	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
383
384	! Ces calculs ne servent pas.
385	! Si necessaire, decommenter ces variables et les calculs...
386	! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
387	! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
388	c
389	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
390	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
391	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
392	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
393
394	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
395	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
396	REAL,ALLOCATABLE,SAVE :: zpk_omp(:,:)
397	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
398	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
399	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
400	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
401	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
402	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
403	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
404	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
405	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
406	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
407	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
408	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
409	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
410
411	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
412	! Introduction du splitting (FH)
413	! Question pour Yann :
414	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
415	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
416	! soit allocatable (plutot par exemple que de passer une dimension
417	! dépendant du process en argument des routines) et que, du coup,
418	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
419	! Tu confirmes ?
420	! J'ai suivi le même principe pour les zdufic_omp
421	! Mais c'est surement bien que tu controles.
422	!
423
424	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
425	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
426	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
427	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
428	REAL jH_cur_split,zdt_split
429	LOGICAL debut_split,lafin_split
430	INTEGER isplit
431	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
432
433	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp,
434	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
435	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
436	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
437	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
438
439	LOGICAL,SAVE :: first_omp=.true.
440	c$OMP THREADPRIVATE(first_omp)
441
442	REAL zsin(iim),zcos(iim),z1(iim)
443	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
444	REAL unskap, pksurcp
445	save unskap
446
447	REAL SSUM
448
449	LOGICAL,SAVE :: firstcal=.true., debut=.true.
450	c$OMP THREADPRIVATE(firstcal,debut)
451
452	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
453	INTEGER :: ierr
454	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
455	INTEGER, dimension(4) :: Req
456	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
457	integer :: k,kstart,kend
458	INTEGER :: offset
459
460	LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
461
462	! For Titan only right now:
463	! to allow for 2D computation of microphys and chemistry
464	LOGICAL,save :: flag_moyzon
465	REAL,allocatable,save :: tmpvar(:,:)
466	REAL,allocatable,save :: tmpvarp1(:,:)
467	REAL,allocatable,save :: tmpvarbar(:)
468	REAL,allocatable,save :: tmpvarbarp1(:)
469	real :: zz1,zz2
470
471	c-----------------------------------------------------------------------
472
473	c 1. Initialisations :
474	c --------------------
475
476
477	klon=klon_mpi
478
479	IF ( firstcal ) THEN
480	debut = .TRUE.
481	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
482	write(lunout,*) 'STOP dans calfis'
483	write(lunout,*)
484	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
485	write(lunout,*) ' ngridmx jjm iim '
486	write(lunout,*) ngridmx,jjm,iim
487	STOP
488	ENDIF
489
490	unskap = 1./ kappa
491
492	c$OMP MASTER
493	flag_moyzon = .false.
494	if(moyzon_ch.or.moyzon_mu) then
495	flag_moyzon = .true.
496	allocate(tmpvar(iip1,llm))
497	allocate(tmpvarp1(iip1,llmp1))
498	allocate(tmpvarbar(llm))
499	allocate(tmpvarbarp1(llmp1))
500	endif
501
502	ALLOCATE(zpsrf(klon))
503	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
504	ALLOCATE(zphi(klon,llm),zphis(klon))
505	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
506	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
507	! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
508	! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
509	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
510	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
511	ALLOCATE(zdpsrf(klon))
512	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
513	ALLOCATE(flxwfi(klon,llm))
514	! ADAPTATION GCM POUR CP(T)
515	ALLOCATE(zteta(klon,llm))
516	ALLOCATE(zpk(klon,llm))
517
518	! zonal means. horizontal dimension should be iip1
519	if (flag_moyzon) call moyzon_init(klon,llm,nqtot)
520
521	c------------------------------------------------------------------
522	c moyennes globales pour les profils de pression et de temperature
523	if(planet_type.eq."titan".or.planet_type.eq."venus") then
524	call AllGather_Field(pp,iip1*jjp1,llmp1)
525	call AllGather_Field(pteta,iip1*jjp1,llm)
526	call AllGather_Field(ppk,iip1*jjp1,llm)
527	call AllGather_Field(pphi,iip1*jjp1,llm)
528	call AllGather_Field(pphis,iip1*jjp1,1)
529	ALLOCATE(plevmoy(llm+1))
530	ALLOCATE(playmoy(llm))
531	ALLOCATE(tmoy(llm))
532	ALLOCATE(tetamoy(llm))
533	ALLOCATE(pkmoy(llm))
534	ALLOCATE(phimoy(0:llm))
535	ALLOCATE(zlevmoy(llm+1))
536	ALLOCATE(zlaymoy(llm))
537	plevmoy=0.
538	do l=1,llmp1
539	do i=1,iip1
540	do j=1,jjp1
541	plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
542	enddo
543	enddo
544	enddo
545	tetamoy=0.
546	pkmoy=0.
547	phimoy=0.
548	do i=1,iip1
549	do j=1,jjp1
550	phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
551	enddo
552	enddo
553	do l=1,llm
554	do i=1,iip1
555	do j=1,jjp1
556	tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
557	pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
558	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
559	enddo
560	enddo
561	enddo
562	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
563	call tpot2t_p(1,llm,tetamoy,tmoy,pkmoy)
564	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
565	zlaymoy(:) = gradrad/(g*rad-phimoy(:))-rad
566	zlevmoy(1) = phimoy(0)/g
567	DO l=2,llm
568	zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
569	zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l))
570	zlevmoy(l)=(zz1zlaymoy(l-1)+zz2zlaymoy(l))/(zz1+zz2)
571	ENDDO
572	zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
573	c-------------------
574	c + lat index
575	allocate(klat(klon))
576	do ig0=1,klon
577	j=index_j(ig0)
578	klat(ig0)=j
579	enddo
580	endif ! planet_type=titan
581	c------------------------------------------------------------------
582
583	c$OMP END MASTER
584	c$OMP BARRIER
585	ELSE
586	debut = .FALSE.
587	ENDIF
588
589
590	c-----------------------------------------------------------------------
591	c 40. transformation des variables dynamiques en variables physiques:
592	c ---------------------------------------------------------------
593
594	c 41. pressions au sol (en Pascals)
595	c ----------------------------------
596
597	c$OMP MASTER
598	call start_timer(timer_physic)
599	c$OMP END MASTER
600
601	c$OMP MASTER
602	!CDIR ON_ADB(index_i)
603	!CDIR ON_ADB(index_j)
604	do ig0=1,klon
605	i=index_i(ig0)
606	j=index_j(ig0)
607	zpsrf(ig0)=pps(i,j)
608	enddo
609	c$OMP END MASTER
610
611
612	c 42. pression intercouches et fonction d'Exner:
613
614	c -----------------------------------------------------------------
615	c .... zplev definis aux (llm +1) interfaces des couches ....
616	c .... zplay definis aux ( llm ) milieux des couches ....
617	c -----------------------------------------------------------------
618
619	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
620
621	c print *,omp_rank,'klon--->',klon
622	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
623	DO l = 1, llmp1
624	!CDIR ON_ADB(index_i)
625	!CDIR ON_ADB(index_j)
626	do ig0=1,klon
627	i=index_i(ig0)
628	j=index_j(ig0)
629	zplev( ig0,l ) = pp(i,j,l)
630	enddo
631	ENDDO
632	c$OMP END DO NOWAIT
633	if (flag_moyzon) then
634	call AllGather_Field(pp,iip1*jjp1,llmp1)
635	j=index_j(1)
636	tmpvarp1(:,:) = pp(:,j,:)
637	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
638	zplevbar_mpi(1,:) = tmpvarbarp1
639	do ig0=2,klon
640	j=index_j(ig0)
641	if (j.ne.index_j(ig0-1)) then
642	tmpvarp1(:,:) = pp(:,j,:)
643	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
644	zplevbar_mpi(ig0,:) = tmpvarbarp1
645	else
646	zplevbar_mpi(ig0,:) = zplevbar_mpi(ig0-1,:)
647	endif
648	enddo
649	endif
650
651	! ADAPTATION GCM POUR CP(T)
652	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
653	DO l=1,llm
654	!CDIR ON_ADB(index_i)
655	!CDIR ON_ADB(index_j)
656	do ig0=1,klon
657	i=index_i(ig0)
658	j=index_j(ig0)
659	zpk(ig0,l)=ppk(i,j,l)
660	zteta(ig0,l)=pteta(i,j,l)
661	enddo
662	ENDDO
663	c$OMP END DO NOWAIT
664	if (flag_moyzon) then
665	call AllGather_Field(pteta,iip1*jjp1,llm)
666	call AllGather_Field(ppk,iip1*jjp1,llm)
667	j=index_j(1)
668	tmpvar(:,:) = pteta(:,j,:)
669	call moyzon(llm,tmpvar,tmpvarbar)
670	ztetabar_mpi(1,:) = tmpvarbar
671	tmpvar(:,:) = ppk(:,j,:)
672	call moyzon(llm,tmpvar,tmpvarbar)
673	zpkbar_mpi(1,:) = tmpvarbar
674	call tpot2t_p(1,llm,ztetabar_mpi(1,:),ztfibar_mpi(1,:),
675	& zpkbar_mpi(1,:))
676	do ig0=2,klon
677	j=index_j(ig0)
678	if (j.ne.index_j(ig0-1)) then
679	tmpvar(:,:) = pteta(:,j,:)
680	call moyzon(llm,tmpvar,tmpvarbar)
681	ztetabar_mpi(ig0,:) = tmpvarbar
682	tmpvar(:,:) = ppk(:,j,:)
683	call moyzon(llm,tmpvar,tmpvarbar)
684	zpkbar_mpi(ig0,:) = tmpvarbar
685	call tpot2t_p(1,llm,ztetabar_mpi(ig0,:),ztfibar_mpi(ig0,:),
686	& zpkbar_mpi(ig0,:))
687	else
688	zpkbar_mpi(ig0,:) = zpkbar_mpi(ig0-1,:)
689	ztetabar_mpi(ig0,:) = ztetabar_mpi(ig0-1,:)
690	ztfibar_mpi(ig0,:) = ztfibar_mpi(ig0-1,:)
691	endif
692	enddo
693	endif
694
695	c 43. temperature naturelle (en K) et pressions milieux couches .
696	c ---------------------------------------------------------------
697
698	! ADAPTATION GCM POUR CP(T)
699	call tpot2t_p(klon,llm,zteta,ztfi,zpk)
700
701	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
702	DO l=1,llm
703	!CDIR ON_ADB(index_i)
704	!CDIR ON_ADB(index_j)
705	do ig0=1,klon
706	i=index_i(ig0)
707	j=index_j(ig0)
708	pksurcp = ppk(i,j,l) / cpp
709	zplay(ig0,l) = preff * pksurcp ** unskap
710	enddo
711	ENDDO
712	c$OMP END DO NOWAIT
713	if (flag_moyzon) then
714	zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap
715	endif
716
717	c 43.bis traceurs (tous intensifs)
718	c ---------------
719
720	DO iq=1,nqtot
721	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
722	DO l=1,llm
723	!CDIR ON_ADB(index_i)
724	!CDIR ON_ADB(index_j)
725	do ig0=1,klon
726	i=index_i(ig0)
727	j=index_j(ig0)
728	zqfi(ig0,l,iq) = pq(i,j,l,iq)
729	enddo
730	ENDDO
731	c$OMP END DO NOWAIT
732	ENDDO ! of DO iq=1,nqtot
733	if (flag_moyzon) then
734	DO iq=1,nqtot
735	call AllGather_Field(pq(:,:,:,iq),iip1*jjp1,llm)
736	j=index_j(1)
737	tmpvar(:,:) = pq(:,j,:,iq)
738	call moyzon(llm,tmpvar,tmpvarbar)
739	zqfibar_mpi(1,:,iq) = tmpvarbar
740	do ig0=2,klon
741	j=index_j(ig0)
742	if (j.ne.index_j(ig0-1)) then
743	tmpvar(:,:) = pq(:,j,:,iq)
744	call moyzon(llm,tmpvar,tmpvarbar)
745	zqfibar_mpi(ig0,:,iq) = tmpvarbar
746	else
747	zqfibar_mpi(ig0,:,iq) = zqfibar_mpi(ig0-1,:,iq)
748	endif
749	enddo
750	ENDDO ! of DO iq=1,nqtot
751	endif
752
753
754	c Geopotentiel calcule par rapport a la surface locale:
755	c -----------------------------------------------------
756
757	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
758
759	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
760
761	c$OMP BARRIER
762	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
763	DO l=1,llm
764	DO ig=1,klon
765	zphi(ig,l)=zphi(ig,l)-zphis(ig)
766	ENDDO
767	ENDDO
768	c$OMP END DO NOWAIT
769
770	if (flag_moyzon) then
771	call AllGather_Field(pphis,iip1*jjp1,1)
772	call AllGather_Field(pphi,iip1*jjp1,llm)
773	j=index_j(1)
774	tmpvar(:,1) = pphis(:,j)
775	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
776	zphisbar_mpi(1) = tmpvarbar(1)
777	tmpvar(:,:) = pphi(:,j,:)
778	call moyzon(llm,tmpvar,tmpvarbar)
779	zphibar_mpi(1,:) = tmpvarbar-zphisbar_mpi(1)
780	do ig0=2,klon
781	j=index_j(ig0)
782	if (j.ne.index_j(ig0-1)) then
783	tmpvar(:,1) = pphis(:,j)
784	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
785	zphisbar_mpi(ig0) = tmpvarbar(1)
786	tmpvar(:,:) = pphi(:,j,:)
787	call moyzon(llm,tmpvar,tmpvarbar)
788	zphibar_mpi(ig0,:) = tmpvarbar-zphisbar_mpi(ig0)
789	else
790	zphisbar_mpi(ig0) = zphisbar_mpi(ig0-1)
791	zphibar_mpi(ig0,:) = zphibar_mpi(ig0-1,:)
792	endif
793	enddo
794	endif
795
796	c
797	c 45. champ u:
798	c ------------
799
800	kstart=1
801	kend=klon
802
803	if (is_north_pole) kstart=2
804	if (is_south_pole) kend=klon-1
805
806	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
807	DO l=1,llm
808	!CDIR ON_ADB(index_i)
809	!CDIR ON_ADB(index_j)
810	!CDIR SPARSE
811	do ig0=kstart,kend
812	i=index_i(ig0)
813	j=index_j(ig0)
814	if (i==1) then
815	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
816	$ + pucov(1,j,l)/cu(1,j) )
817	else
818	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
819	$ + pucov(i,j,l)/cu(i,j) )
820	endif
821	enddo
822	ENDDO
823	c$OMP END DO NOWAIT
824
825	c 46.champ v:
826	c -----------
827
828	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
829	DO l=1,llm
830	!CDIR ON_ADB(index_i)
831	!CDIR ON_ADB(index_j)
832	DO ig0=kstart,kend
833	i=index_i(ig0)
834	j=index_j(ig0)
835	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
836	$ + pvcov(i,j,l)/cv(i,j) )
837
838	ENDDO
839	ENDDO
840	c$OMP END DO NOWAIT
841
842	c 47. champs de vents aux pole nord
843	c ------------------------------
844	c U = 1 / pi * integrale [ v * cos(long) * d long ]
845	c V = 1 / pi * integrale [ v * sin(long) * d long ]
846
847	if (is_north_pole) then
848	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
849	DO l=1,llm
850
851	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
852	DO i=2,iim
853	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
854	ENDDO
855
856	DO i=1,iim
857	zcos(i) = COS(rlonv(i))*z1(i)
858	zsin(i) = SIN(rlonv(i))*z1(i)
859	ENDDO
860
861	zufi(1,l) = SSUM(iim,zcos,1)/pi
862	zvfi(1,l) = SSUM(iim,zsin,1)/pi
863
864	ENDDO
865	c$OMP END DO NOWAIT
866	endif
867
868
869	c 48. champs de vents aux pole sud:
870	c ---------------------------------
871	c U = 1 / pi * integrale [ v * cos(long) * d long ]
872	c V = 1 / pi * integrale [ v * sin(long) * d long ]
873
874	if (is_south_pole) then
875	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
876	DO l=1,llm
877
878	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
879	DO i=2,iim
880	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
881	ENDDO
882
883	DO i=1,iim
884	zcos(i) = COS(rlonv(i))*z1(i)
885	zsin(i) = SIN(rlonv(i))*z1(i)
886	ENDDO
887
888	zufi(klon,l) = SSUM(iim,zcos,1)/pi
889	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
890	ENDDO
891	c$OMP END DO NOWAIT
892	endif
893
894	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
895	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
896
897	c-----------------------------------------------------------------------
898	c Appel de la physique:
899	c ---------------------
900
901	! Appel de la physique: pose probleme quand on tourne
902	! SANS physique, car physiq.F est dans le repertoire phy[]...
903	! Il faut une cle 1
904
905	! Le fait que les arguments de physiq soient differents selon les planetes
906	! ne pose pas de probleme a priori.
907
908
909	c$OMP BARRIER
910	if (first_omp) then
911	klon=klon_omp
912
913	allocate(zplev_omp(klon,llm+1))
914	allocate(zplay_omp(klon,llm))
915	allocate(zpk_omp(klon,llm))
916	allocate(zphi_omp(klon,llm))
917	allocate(zphis_omp(klon))
918	allocate(presnivs_omp(llm))
919	allocate(zufi_omp(klon,llm))
920	allocate(zvfi_omp(klon,llm))
921	allocate(ztfi_omp(klon,llm))
922	allocate(zqfi_omp(klon,llm,nqtot))
923	allocate(zdufi_omp(klon,llm))
924	allocate(zdvfi_omp(klon,llm))
925	allocate(zdtfi_omp(klon,llm))
926	allocate(zdqfi_omp(klon,llm,nqtot))
927	allocate(zdufic_omp(klon,llm))
928	allocate(zdvfic_omp(klon,llm))
929	allocate(zdtfic_omp(klon,llm))
930	allocate(zdqfic_omp(klon,llm,nqtot))
931	allocate(zdpsrf_omp(klon))
932	allocate(flxwfi_omp(klon,llm))
933
934	if (flag_moyzon) call moyzon_init_omp(klon,llm,nqtot)
935
936	first_omp=.false.
937	endif
938
939
940	klon=klon_omp
941	offset=klon_omp_begin-1
942
943	do l=1,llm+1
944	do i=1,klon
945	zplev_omp(i,l)=zplev(offset+i,l)
946	enddo
947	enddo
948
949	do l=1,llm
950	do i=1,klon
951	zplay_omp(i,l)=zplay(offset+i,l)
952	enddo
953	enddo
954
955	do l=1,llm
956	do i=1,klon
957	zpk_omp(i,l)=zpk(offset+i,l)
958	enddo
959	enddo
960
961	do l=1,llm
962	do i=1,klon
963	zphi_omp(i,l)=zphi(offset+i,l)
964	enddo
965	enddo
966
967	do i=1,klon
968	zphis_omp(i)=zphis(offset+i)
969	enddo
970
971
972	do l=1,llm
973	presnivs_omp(l)=presnivs(l)
974	enddo
975
976	do l=1,llm
977	do i=1,klon
978	zufi_omp(i,l)=zufi(offset+i,l)
979	enddo
980	enddo
981
982	do l=1,llm
983	do i=1,klon
984	zvfi_omp(i,l)=zvfi(offset+i,l)
985	enddo
986	enddo
987
988	do l=1,llm
989	do i=1,klon
990	ztfi_omp(i,l)=ztfi(offset+i,l)
991	enddo
992	enddo
993
994	do iq=1,nqtot
995	do l=1,llm
996	do i=1,klon
997	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
998	enddo
999	enddo
1000	enddo
1001
1002	do l=1,llm
1003	do i=1,klon
1004	zdufi_omp(i,l)=zdufi(offset+i,l)
1005	enddo
1006	enddo
1007
1008	do l=1,llm
1009	do i=1,klon
1010	zdvfi_omp(i,l)=zdvfi(offset+i,l)
1011	enddo
1012	enddo
1013
1014	do l=1,llm
1015	do i=1,klon
1016	zdtfi_omp(i,l)=zdtfi(offset+i,l)
1017	enddo
1018	enddo
1019
1020	do iq=1,nqtot
1021	do l=1,llm
1022	do i=1,klon
1023	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
1024	enddo
1025	enddo
1026	enddo
1027
1028	do i=1,klon
1029	zdpsrf_omp(i)=zdpsrf(offset+i)
1030	enddo
1031
1032	do l=1,llm
1033	do i=1,klon
1034	flxwfi_omp(i,l)=flxwfi(offset+i,l)
1035	enddo
1036	enddo
1037
1038	if (flag_moyzon) then
1039	do l=1,llm+1
1040	do i=1,klon
1041	zplevbar(i,l)=zplevbar_mpi(offset+i,l)
1042	enddo
1043	enddo
1044
1045	do l=1,llm
1046	do i=1,klon
1047	zplaybar(i,l)=zplaybar_mpi(offset+i,l)
1048	enddo
1049	enddo
1050
1051	do l=1,llm
1052	do i=1,klon
1053	zphibar(i,l)=zphibar_mpi(offset+i,l)
1054	enddo
1055	enddo
1056
1057	do i=1,klon
1058	zphisbar(i)=zphisbar_mpi(offset+i)
1059	enddo
1060
1061	do l=1,llm
1062	do i=1,klon
1063	ztfibar(i,l)=ztfibar_mpi(offset+i,l)
1064	enddo
1065	enddo
1066
1067	do iq=1,nqtot
1068	do l=1,llm
1069	do i=1,klon
1070	zqfibar(i,l,iq)=zqfibar_mpi(offset+i,l,iq)
1071	enddo
1072	enddo
1073	enddo
1074	endif
1075
1076
1077	c$OMP BARRIER
1078
1079	!$OMP MASTER
1080	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
1081	!$OMP END MASTER
1082	zdt_split=dtphys/nsplit_phys
1083	zdufic_omp(:,:)=0.
1084	zdvfic_omp(:,:)=0.
1085	zdtfic_omp(:,:)=0.
1086	zdqfic_omp(:,:,:)=0.
1087
1088	do isplit=1,nsplit_phys
1089
1090	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
1091	debut_split=debut.and.isplit==1
1092	lafin_split=lafin.and.isplit==nsplit_phys
1093
1094
1095	if (planet_type=="earth") then
1096	CALL physiq (klon,
1097	. llm,
1098	. debut_split,
1099	. lafin_split,
1100	. jD_cur,
1101	. jH_cur_split,
1102	. zdt_split,
1103	. zplev_omp,
1104	. zplay_omp,
1105	. zphi_omp,
1106	. zphis_omp,
1107	. presnivs_omp,
1108	. zufi_omp,
1109	. zvfi_omp,
1110	. ztfi_omp,
1111	. zqfi_omp,
1112	c#ifdef INCA
1113	. flxwfi_omp,
1114	c#endif
1115	. zdufi_omp,
1116	. zdvfi_omp,
1117	. zdtfi_omp,
1118	. zdqfi_omp,
1119	. zdpsrf_omp,
1120	. pducov)
1121
1122	else if ( planet_type=="generic" ) then
1123
1124	CALL physiq (klon, !! ngrid
1125	. llm, !! nlayer
1126	. nqtot, !! nq
1127	. tname, !! tracer names from dynamical core (given in infotrac)
1128	. debut_split, !! firstcall
1129	. lafin_split, !! lastcall
1130	. jD_cur, !! pday. see leapfrog_p
1131	. jH_cur_split, !! ptime "fraction of day"
1132	. zdt_split, !! ptimestep
1133	. zplev_omp, !! pplev
1134	. zplay_omp, !! pplay
1135	. zphi_omp, !! pphi
1136	. zufi_omp, !! pu
1137	. zvfi_omp, !! pv
1138	. ztfi_omp, !! pt
1139	. zqfi_omp, !! pq
1140	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
1141	. zdufi_omp, !! pdu
1142	. zdvfi_omp, !! pdv
1143	. zdtfi_omp, !! pdt
1144	. zdqfi_omp, !! pdq
1145	. zdpsrf_omp, !! pdpsrf
1146	. tracerdyn) !! tracerdyn <-- utilite ???
1147
1148	else if ( planet_type=="mars" ) then
1149
1150	CALL physiq (klon, ! ngrid
1151	. llm, ! nlayer
1152	. nqtot, ! nq
1153	. debut_split, ! firstcall
1154	. lafin_split, ! lastcall
1155	. jD_cur, ! pday
1156	. jH_cur_split, ! ptime
1157	. zdt_split, ! ptimestep
1158	. zplev_omp, ! pplev
1159	. zplay_omp, ! pplay
1160	. zphi_omp, ! pphi
1161	. zufi_omp, ! pu
1162	. zvfi_omp, ! pv
1163	. ztfi_omp, ! pt
1164	. zqfi_omp, ! pq
1165	. flxwfi_omp, ! pw
1166	. zdufi_omp, ! pdu
1167	. zdvfi_omp, ! pdv
1168	. zdtfi_omp, ! pdt
1169	. zdqfi_omp, ! pdq
1170	. zdpsrf_omp, ! pdpsrf
1171	. tracerdyn) ! tracerdyn (somewhat obsolete)
1172
1173	else if ((planet_type=="titan").or.(planet_type=="venus")) then
1174
1175	CALL physiq (klon,
1176	. llm,
1177	. nqtot,
1178	. debut_split,
1179	. lafin_split,
1180	. jD_cur,
1181	. jH_cur_split,
1182	. zdt_split,
1183	. zplev_omp,
1184	. zplay_omp,
1185	. zpk_omp,
1186	. zphi_omp,
1187	. zphis_omp,
1188	. presnivs_omp,
1189	. zufi_omp,
1190	. zvfi_omp,
1191	. ztfi_omp,
1192	. zqfi_omp,
1193	. flxwfi_omp,
1194	. zdufi_omp,
1195	. zdvfi_omp,
1196	. zdtfi_omp,
1197	. zdqfi_omp,
1198	. zdpsrf_omp)
1199
1200	else ! unknown "planet_type"
1201
1202	write(lunout,*) "calfis_p: error, unknown planet_type: ",
1203	& trim(planet_type)
1204	stop
1205
1206	endif ! planet_type
1207	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
1208	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
1209	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
1210	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
1211
1212	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
1213	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
1214	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
1215	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
1216
1217	enddo
1218
1219	! of #ifdef 1
1220
1221	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
1222	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
1223	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
1224	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
1225
1226	c$OMP BARRIER
1227
1228	do l=1,llm+1
1229	do i=1,klon
1230	zplev(offset+i,l)=zplev_omp(i,l)
1231	enddo
1232	enddo
1233
1234	do l=1,llm
1235	do i=1,klon
1236	zplay(offset+i,l)=zplay_omp(i,l)
1237	enddo
1238	enddo
1239
1240	do l=1,llm
1241	do i=1,klon
1242	zphi(offset+i,l)=zphi_omp(i,l)
1243	enddo
1244	enddo
1245
1246
1247	do i=1,klon
1248	zphis(offset+i)=zphis_omp(i)
1249	enddo
1250
1251
1252	do l=1,llm
1253	presnivs(l)=presnivs_omp(l)
1254	enddo
1255
1256	do l=1,llm
1257	do i=1,klon
1258	zufi(offset+i,l)=zufi_omp(i,l)
1259	enddo
1260	enddo
1261
1262	do l=1,llm
1263	do i=1,klon
1264	zvfi(offset+i,l)=zvfi_omp(i,l)
1265	enddo
1266	enddo
1267
1268	do l=1,llm
1269	do i=1,klon
1270	ztfi(offset+i,l)=ztfi_omp(i,l)
1271	enddo
1272	enddo
1273
1274	do iq=1,nqtot
1275	do l=1,llm
1276	do i=1,klon
1277	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
1278	enddo
1279	enddo
1280	enddo
1281
1282	do l=1,llm
1283	do i=1,klon
1284	zdufi(offset+i,l)=zdufi_omp(i,l)
1285	enddo
1286	enddo
1287
1288	do l=1,llm
1289	do i=1,klon
1290	zdvfi(offset+i,l)=zdvfi_omp(i,l)
1291	enddo
1292	enddo
1293
1294	do l=1,llm
1295	do i=1,klon
1296	zdtfi(offset+i,l)=zdtfi_omp(i,l)
1297	enddo
1298	enddo
1299
1300	do iq=1,nqtot
1301	do l=1,llm
1302	do i=1,klon
1303	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
1304	enddo
1305	enddo
1306	enddo
1307
1308	do i=1,klon
1309	zdpsrf(offset+i)=zdpsrf_omp(i)
1310	enddo
1311
1312
1313	klon=klon_mpi
1314	500 CONTINUE
1315	c$OMP BARRIER
1316
1317	c$OMP MASTER
1318	call stop_timer(timer_physic)
1319	c$OMP END MASTER
1320
1321	IF (using_mpi) THEN
1322
1323	if (MPI_rank>0) then
1324
1325	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1326	DO l=1,llm
1327	du_send(1:iim,l)=zdufi(1:iim,l)
1328	dv_send(1:iim,l)=zdvfi(1:iim,l)
1329	ENDDO
1330	c$OMP END DO NOWAIT
1331
1332	c$OMP BARRIER
1333	c$OMP MASTER
1334	!$OMP CRITICAL (MPI)
1335	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
1336	& COMM_LMDZ,Req(1),ierr)
1337	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
1338	& COMM_LMDZ,Req(2),ierr)
1339	!$OMP END CRITICAL (MPI)
1340	c$OMP END MASTER
1341	c$OMP BARRIER
1342
1343	endif
1344
1345	if (MPI_rank<MPI_Size-1) then
1346	c$OMP BARRIER
1347	c$OMP MASTER
1348	!$OMP CRITICAL (MPI)
1349	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
1350	& COMM_LMDZ,Req(3),ierr)
1351	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
1352	& COMM_LMDZ,Req(4),ierr)
1353	!$OMP END CRITICAL (MPI)
1354	c$OMP END MASTER
1355	endif
1356
1357	c$OMP BARRIER
1358
1359
1360	c$OMP MASTER
1361	!$OMP CRITICAL (MPI)
1362	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
1363	call MPI_WAITALL(4,Req(1),Status,ierr)
1364	else if (MPI_rank>0) then
1365	call MPI_WAITALL(2,Req(1),Status,ierr)
1366	else if (MPI_rank <MPI_Size-1) then
1367	call MPI_WAITALL(2,Req(3),Status,ierr)
1368	endif
1369	!$OMP END CRITICAL (MPI)
1370	c$OMP END MASTER
1371
1372	c$OMP BARRIER
1373
1374	ENDIF ! using_mpi
1375
1376
1377	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1378	DO l=1,llm
1379
1380	zdufi2(1:klon,l)=zdufi(1:klon,l)
1381	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
1382
1383	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
1384	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
1385
1386	pdhfi(:,jj_begin,l)=0
1387	pdqfi(:,jj_begin,l,:)=0
1388	pdufi(:,jj_begin,l)=0
1389	pdvfi(:,jj_begin,l)=0
1390
1391	if (.not. is_south_pole) then
1392	pdhfi(:,jj_end,l)=0
1393	pdqfi(:,jj_end,l,:)=0
1394	pdufi(:,jj_end,l)=0
1395	pdvfi(:,jj_end,l)=0
1396	endif
1397
1398	ENDDO
1399	c$OMP END DO NOWAIT
1400
1401	c$OMP MASTER
1402	pdpsfi(:,jj_begin)=0
1403	if (.not. is_south_pole) then
1404	pdpsfi(:,jj_end)=0
1405	endif
1406	c$OMP END MASTER
1407	c-----------------------------------------------------------------------
1408	c transformation des tendances physiques en tendances dynamiques:
1409	c ---------------------------------------------------------------
1410
1411	c tendance sur la pression :
1412	c -----------------------------------
1413	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
1414	c
1415	c 62. enthalpie potentielle
1416	c ---------------------
1417
1418	kstart=1
1419	kend=klon
1420
1421	if (is_north_pole) kstart=2
1422	if (is_south_pole) kend=klon-1
1423
1424	! ADAPTATION GCM POUR CP(T)
1425	call t2tpot_p(klon,llm,ztfi,zteta,zpk)
1426
1427
1428	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1429	DO l=1,llm
1430	!CDIR ON_ADB(index_i)
1431	!CDIR ON_ADB(index_j)
1432	!cdir NODEP
1433	do ig0=kstart,kend
1434	i=index_i(ig0)
1435	j=index_j(ig0)
1436	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1437	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1438	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1439	if (i==1) then
1440	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1441	endif
1442	enddo
1443
1444	if (is_north_pole) then
1445	DO i=1,iip1
1446	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
1447	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
1448	enddo
1449	endif
1450
1451	if (is_south_pole) then
1452	DO i=1,iip1
1453	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1454	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
1455	ENDDO
1456	endif
1457	ENDDO
1458	c$OMP END DO NOWAIT
1459
1460	c 62. humidite specifique
1461	c ---------------------
1462	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1463	! DO iq=1,nqtot
1464	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1465	! DO l=1,llm
1466	!!!cdir NODEP
1467	! do ig0=kstart,kend
1468	! i=index_i(ig0)
1469	! j=index_j(ig0)
1470	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1471	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1472	! enddo
1473	!
1474	! if (is_north_pole) then
1475	! do i=1,iip1
1476	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1477	! enddo
1478	! endif
1479	!
1480	! if (is_south_pole) then
1481	! do i=1,iip1
1482	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1483	! enddo
1484	! endif
1485	! ENDDO
1486	!c$OMP END DO NOWAIT
1487	! ENDDO
1488
1489	c 63. traceurs (tous en intensifs)
1490	c ------------
1491	C initialisation des tendances
1492
1493	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1494	DO l=1,llm
1495	pdqfi(:,:,l,:)=0.
1496	ENDDO
1497	c$OMP END DO NOWAIT
1498
1499	C
1500	!cdir NODEP
1501	DO iq=1,nqtot
1502	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1503	DO l=1,llm
1504	!CDIR ON_ADB(index_i)
1505	!CDIR ON_ADB(index_j)
1506	!cdir NODEP
1507	DO ig0=kstart,kend
1508	i=index_i(ig0)
1509	j=index_j(ig0)
1510	pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1511	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1512	ENDDO
1513
1514	IF (is_north_pole) then
1515	DO i=1,iip1
1516	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1517	ENDDO
1518	ENDIF
1519
1520	IF (is_south_pole) then
1521	DO i=1,iip1
1522	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1523	ENDDO
1524	ENDIF
1525
1526	ENDDO
1527	c$OMP END DO NOWAIT
1528	ENDDO
1529
1530	c 65. champ u:
1531	c ------------
1532	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1533	DO l=1,llm
1534	!CDIR ON_ADB(index_i)
1535	!CDIR ON_ADB(index_j)
1536	!cdir NODEP
1537	do ig0=kstart,kend
1538	i=index_i(ig0)
1539	j=index_j(ig0)
1540
1541	if (i/=iim) then
1542	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1543	endif
1544
1545	if (i==1) then
1546	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1547	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1548	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1549	endif
1550
1551	enddo
1552
1553	if (is_north_pole) then
1554	DO i=1,iip1
1555	pdufi(i,1,l) = 0.
1556	ENDDO
1557	endif
1558
1559	if (is_south_pole) then
1560	DO i=1,iip1
1561	pdufi(i,jjp1,l) = 0.
1562	ENDDO
1563	endif
1564
1565	ENDDO
1566	c$OMP END DO NOWAIT
1567
1568	c 67. champ v:
1569	c ------------
1570
1571	kstart=1
1572	kend=klon
1573
1574	if (is_north_pole) kstart=2
1575	if (is_south_pole) kend=klon-1-iim
1576
1577	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1578	DO l=1,llm
1579	!CDIR ON_ADB(index_i)
1580	!CDIR ON_ADB(index_j)
1581	!cdir NODEP
1582	do ig0=kstart,kend
1583	i=index_i(ig0)
1584	j=index_j(ig0)
1585	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1586	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1587	$ zdvfi2(ig0+iim,l))
1588	$ *cv(i,j)
1589	enddo
1590
1591	ENDDO
1592	c$OMP END DO NOWAIT
1593
1594
1595	c 68. champ v pres des poles:
1596	c ---------------------------
1597	c v = U * cos(long) + V * SIN(long)
1598
1599	if (is_north_pole) then
1600
1601	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1602	DO l=1,llm
1603
1604	DO i=1,iim
1605	pdvfi(i,1,l)=
1606	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1607
1608	pdvfi(i,1,l)=
1609	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1610	ENDDO
1611
1612	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1613
1614	ENDDO
1615	c$OMP END DO NOWAIT
1616
1617	endif
1618
1619	if (is_south_pole) then
1620
1621	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1622	DO l=1,llm
1623
1624	DO i=1,iim
1625	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1626	$ +zdvfi(klon,l)*SIN(rlonv(i))
1627
1628	pdvfi(i,jjm,l)=
1629	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1630	ENDDO
1631
1632	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1633
1634	ENDDO
1635	c$OMP END DO NOWAIT
1636
1637	endif
1638	c-----------------------------------------------------------------------
1639
1640	700 CONTINUE
1641
1642	firstcal = .FALSE.
1643
1644	! of #ifdef 1
1645	RETURN
1646	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats: