Context Navigation

calfis_p.F @ 1017

Last change on this file since 1017 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

File size: 31.6 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	#ifdef CPP_PHYS
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	#endif
35	USE parallel, ONLY : omp_chunk, using_mpi
36	USE Write_Field
37	Use Write_field_p
38	USE Times
39	USE infotrac
40	USE control_mod
41	USE cpdet_mod, only: tpot2t_p, t2tpot_p
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 (K/s)
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	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
129	! qui lui meme ne sert a rien dans la routine telle qu'elle est
130	! ecrite, et que j'ai donc commente....
131	REAL pdq(iip1,jjp1,llm,nqtot)
132	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
133	c
134	REAL pps(iip1,jjp1)
135	REAL pp(iip1,jjp1,llmp1)
136	REAL ppk(iip1,jjp1,llm)
137	c
138	REAL pdvfi(iip1,jjm,llm)
139	REAL pdufi(iip1,jjp1,llm)
140	REAL pdhfi(iip1,jjp1,llm)
141	REAL pdqfi(iip1,jjp1,llm,nqtot)
142	REAL pdpsfi(iip1,jjp1)
143
144	#ifdef CPP_PHYS
145	c Local variables :
146	c -----------------
147
148	INTEGER i,j,l,ig0,ig,iq,iiq
149	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
150	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
151	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
152	c
153	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
154	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
155	! ADAPTATION GCM POUR CP(T)
156	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
157	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
158	c
159	! Ces calculs ne servent pas.
160	! Si necessaire, decommenter ces variables et les calculs...
161	! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
162	! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
163	c
164	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
165	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
166	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
167	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
168
169	c
170	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
171	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
172	REAL,ALLOCATABLE,SAVE :: zpk_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,zpk_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 ! for generic/mars physics call ; possibly to get rid of
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	! ADAPTATION GCM POUR CP(T)
286	ALLOCATE(zteta(klon,llm))
287	ALLOCATE(zpk(klon,llm))
288	c$OMP END MASTER
289	c$OMP BARRIER
290	ELSE
291	debut = .FALSE.
292	ENDIF
293
294	c
295	c
296	c-----------------------------------------------------------------------
297	c 40. transformation des variables dynamiques en variables physiques:
298	c ---------------------------------------------------------------
299
300	c 41. pressions au sol (en Pascals)
301	c ----------------------------------
302
303	c$OMP MASTER
304	call start_timer(timer_physic)
305	c$OMP END MASTER
306
307	c$OMP MASTER
308	!CDIR ON_ADB(index_i)
309	!CDIR ON_ADB(index_j)
310	do ig0=1,klon
311	i=index_i(ig0)
312	j=index_j(ig0)
313	zpsrf(ig0)=pps(i,j)
314	enddo
315	c$OMP END MASTER
316
317
318	c 42. pression intercouches et fonction d'Exner:
319	c
320	c -----------------------------------------------------------------
321	c .... zplev definis aux (llm +1) interfaces des couches ....
322	c .... zplay definis aux ( llm ) milieux des couches ....
323	c -----------------------------------------------------------------
324
325	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
326	c
327	unskap = 1./ kappa
328	c
329	c print *,omp_rank,'klon--->',klon
330	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
331	DO l = 1, llmp1
332	!CDIR ON_ADB(index_i)
333	!CDIR ON_ADB(index_j)
334	do ig0=1,klon
335	i=index_i(ig0)
336	j=index_j(ig0)
337	zplev( ig0,l ) = pp(i,j,l)
338	enddo
339	ENDDO
340	c$OMP END DO NOWAIT
341	! ADAPTATION GCM POUR CP(T)
342	DO l=1,llm
343	!CDIR ON_ADB(index_i)
344	!CDIR ON_ADB(index_j)
345	do ig0=1,klon
346	i=index_i(ig0)
347	j=index_j(ig0)
348	zpk(ig0,l)=ppk(i,j,l)
349	zteta(ig0,l)=pteta(i,j,l)
350	enddo
351	ENDDO
352	c$OMP END DO NOWAIT
353
354	c
355	c
356
357	c 43. temperature naturelle (en K) et pressions milieux couches .
358	c ---------------------------------------------------------------
359
360	! ADAPTATION GCM POUR CP(T)
361	call tpot2t_p(klon,llm,zteta,ztfi,zpk)
362
363	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
364	DO l=1,llm
365	!CDIR ON_ADB(index_i)
366	!CDIR ON_ADB(index_j)
367	do ig0=1,klon
368	i=index_i(ig0)
369	j=index_j(ig0)
370	pksurcp = ppk(i,j,l) / cpp
371	zplay(ig0,l) = preff * pksurcp ** unskap
372	! ztfi(ig0,l) = pteta(i,j,l) * pksurcp
373	enddo
374
375	ENDDO
376	c$OMP END DO NOWAIT
377
378	c 43.bis traceurs (tous intensifs)
379	c ---------------
380	c
381
382	DO iq=1,nqtot
383	iiq=niadv(iq)
384	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
385	DO l=1,llm
386	!CDIR ON_ADB(index_i)
387	!CDIR ON_ADB(index_j)
388	do ig0=1,klon
389	i=index_i(ig0)
390	j=index_j(ig0)
391	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
392	enddo
393	ENDDO
394	c$OMP END DO NOWAIT
395	ENDDO ! of DO iq=1,nqtot
396
397
398	c Geopotentiel calcule par rapport a la surface locale:
399	c -----------------------------------------------------
400
401	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
402
403	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
404
405	c$OMP BARRIER
406
407	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
408	DO l=1,llm
409	DO ig=1,klon
410	zphi(ig,l)=zphi(ig,l)-zphis(ig)
411	ENDDO
412	ENDDO
413	c$OMP END DO NOWAIT
414
415
416	c
417	c 45. champ u:
418	c ------------
419
420	kstart=1
421	kend=klon
422
423	if (is_north_pole) kstart=2
424	if (is_south_pole) kend=klon-1
425
426	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
427	DO l=1,llm
428	!CDIR ON_ADB(index_i)
429	!CDIR ON_ADB(index_j)
430	!CDIR SPARSE
431	do ig0=kstart,kend
432	i=index_i(ig0)
433	j=index_j(ig0)
434	if (i==1) then
435	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
436	$ + pucov(1,j,l)/cu(1,j) )
437	else
438	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
439	$ + pucov(i,j,l)/cu(i,j) )
440	endif
441	enddo
442	ENDDO
443	c$OMP END DO NOWAIT
444
445	c 46.champ v:
446	c -----------
447
448	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
449	DO l=1,llm
450	!CDIR ON_ADB(index_i)
451	!CDIR ON_ADB(index_j)
452	DO ig0=kstart,kend
453	i=index_i(ig0)
454	j=index_j(ig0)
455	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
456	$ + pvcov(i,j,l)/cv(i,j) )
457
458	ENDDO
459	ENDDO
460	c$OMP END DO NOWAIT
461
462	c 47. champs de vents aux pole nord
463	c ------------------------------
464	c U = 1 / pi * integrale [ v * cos(long) * d long ]
465	c V = 1 / pi * integrale [ v * sin(long) * d long ]
466
467	if (is_north_pole) then
468	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
469	DO l=1,llm
470
471	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
472	DO i=2,iim
473	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
474	ENDDO
475
476	DO i=1,iim
477	zcos(i) = COS(rlonv(i))*z1(i)
478	zsin(i) = SIN(rlonv(i))*z1(i)
479	ENDDO
480
481	zufi(1,l) = SSUM(iim,zcos,1)/pi
482	zvfi(1,l) = SSUM(iim,zsin,1)/pi
483
484	ENDDO
485	c$OMP END DO NOWAIT
486	endif
487
488
489	c 48. champs de vents aux pole sud:
490	c ---------------------------------
491	c U = 1 / pi * integrale [ v * cos(long) * d long ]
492	c V = 1 / pi * integrale [ v * sin(long) * d long ]
493
494	if (is_south_pole) then
495	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
496	DO l=1,llm
497
498	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
499	DO i=2,iim
500	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
501	ENDDO
502
503	DO i=1,iim
504	zcos(i) = COS(rlonv(i))*z1(i)
505	zsin(i) = SIN(rlonv(i))*z1(i)
506	ENDDO
507
508	zufi(klon,l) = SSUM(iim,zcos,1)/pi
509	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
510	ENDDO
511	c$OMP END DO NOWAIT
512	endif
513
514
515	IF (is_sequential.and.(planet_type=="earth")) THEN
516	#ifdef CPP_EARTH
517	! PVtheta calls tetalevel, which is in the physics
518	cIM calcul PV a teta=350, 380, 405K
519	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
520	$ ztfi,zplay,zplev,
521	$ ntetaSTD,rtetaSTD,PVteta)
522	c
523	#endif
524	ENDIF
525
526	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
527	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
528
529	c-----------------------------------------------------------------------
530	c Appel de la physique:
531	c ---------------------
532
533	! Appel de la physique: pose probleme quand on tourne
534	! SANS physique, car physiq.F est dans le repertoire phy[]...
535	! Il faut une cle CPP_PHYS
536
537	! Le fait que les arguments de physiq soient differents selon les planetes
538	! ne pose pas de probleme a priori.
539
540
541	c$OMP BARRIER
542	if (first_omp) then
543	klon=klon_omp
544
545	allocate(zplev_omp(klon,llm+1))
546	allocate(zplay_omp(klon,llm))
547	allocate(zpk_omp(klon,llm))
548	allocate(zphi_omp(klon,llm))
549	allocate(zphis_omp(klon))
550	allocate(presnivs_omp(llm))
551	allocate(zufi_omp(klon,llm))
552	allocate(zvfi_omp(klon,llm))
553	allocate(ztfi_omp(klon,llm))
554	allocate(zqfi_omp(klon,llm,nqtot))
555	allocate(zdufi_omp(klon,llm))
556	allocate(zdvfi_omp(klon,llm))
557	allocate(zdtfi_omp(klon,llm))
558	allocate(zdqfi_omp(klon,llm,nqtot))
559	allocate(zdufic_omp(klon,llm))
560	allocate(zdvfic_omp(klon,llm))
561	allocate(zdtfic_omp(klon,llm))
562	allocate(zdqfic_omp(klon,llm,nqtot))
563	allocate(zdpsrf_omp(klon))
564	allocate(flxwfi_omp(klon,llm))
565	first_omp=.false.
566	endif
567
568
569	klon=klon_omp
570	offset=klon_omp_begin-1
571
572	do l=1,llm+1
573	do i=1,klon
574	zplev_omp(i,l)=zplev(offset+i,l)
575	enddo
576	enddo
577
578	do l=1,llm
579	do i=1,klon
580	zplay_omp(i,l)=zplay(offset+i,l)
581	enddo
582	enddo
583
584	do l=1,llm
585	do i=1,klon
586	zpk_omp(i,l)=zpk(offset+i,l)
587	enddo
588	enddo
589
590	do l=1,llm
591	do i=1,klon
592	zphi_omp(i,l)=zphi(offset+i,l)
593	enddo
594	enddo
595
596	do i=1,klon
597	zphis_omp(i)=zphis(offset+i)
598	enddo
599
600
601	do l=1,llm
602	presnivs_omp(l)=presnivs(l)
603	enddo
604
605	do l=1,llm
606	do i=1,klon
607	zufi_omp(i,l)=zufi(offset+i,l)
608	enddo
609	enddo
610
611	do l=1,llm
612	do i=1,klon
613	zvfi_omp(i,l)=zvfi(offset+i,l)
614	enddo
615	enddo
616
617	do l=1,llm
618	do i=1,klon
619	ztfi_omp(i,l)=ztfi(offset+i,l)
620	enddo
621	enddo
622
623	do iq=1,nqtot
624	do l=1,llm
625	do i=1,klon
626	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
627	enddo
628	enddo
629	enddo
630
631	do l=1,llm
632	do i=1,klon
633	zdufi_omp(i,l)=zdufi(offset+i,l)
634	enddo
635	enddo
636
637	do l=1,llm
638	do i=1,klon
639	zdvfi_omp(i,l)=zdvfi(offset+i,l)
640	enddo
641	enddo
642
643	do l=1,llm
644	do i=1,klon
645	zdtfi_omp(i,l)=zdtfi(offset+i,l)
646	enddo
647	enddo
648
649	do iq=1,nqtot
650	do l=1,llm
651	do i=1,klon
652	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
653	enddo
654	enddo
655	enddo
656
657	do i=1,klon
658	zdpsrf_omp(i)=zdpsrf(offset+i)
659	enddo
660
661	do l=1,llm
662	do i=1,klon
663	flxwfi_omp(i,l)=flxwfi(offset+i,l)
664	enddo
665	enddo
666
667	c$OMP BARRIER
668
669	!$OMP MASTER
670	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
671	!$OMP END MASTER
672	zdt_split=dtphys/nsplit_phys
673	zdufic_omp(:,:)=0.
674	zdvfic_omp(:,:)=0.
675	zdtfic_omp(:,:)=0.
676	zdqfic_omp(:,:,:)=0.
677
678	#ifdef CPP_PHYS
679	do isplit=1,nsplit_phys
680
681	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
682	debut_split=debut.and.isplit==1
683	lafin_split=lafin.and.isplit==nsplit_phys
684
685
686	if (planet_type=="earth") then
687	CALL physiq (klon,
688	. llm,
689	. debut_split,
690	. lafin_split,
691	. jD_cur,
692	. jH_cur_split,
693	. zdt_split,
694	. zplev_omp,
695	. zplay_omp,
696	. zphi_omp,
697	. zphis_omp,
698	. presnivs_omp,
699	. zufi_omp,
700	. zvfi_omp,
701	. ztfi_omp,
702	. zqfi_omp,
703	c#ifdef INCA
704	. flxwfi_omp,
705	c#endif
706	. zdufi_omp,
707	. zdvfi_omp,
708	. zdtfi_omp,
709	. zdqfi_omp,
710	. zdpsrf_omp,
711	cIM diagnostique PVteta, Amip2
712	. pducov,
713	. PVteta)
714
715	else if ( planet_type=="generic" ) then
716
717	CALL physiq (klon, !! ngrid
718	. llm, !! nlayer
719	. nqtot, !! nq
720	. tname, !! tracer names from dynamical core (given in infotrac)
721	. debut_split, !! firstcall
722	. lafin_split, !! lastcall
723	. jD_cur, !! pday. see leapfrog_p
724	. jH_cur_split, !! ptime "fraction of day"
725	. zdt_split, !! ptimestep
726	. zplev_omp, !! pplev
727	. zplay_omp, !! pplay
728	. zphi_omp, !! pphi
729	. zufi_omp, !! pu
730	. zvfi_omp, !! pv
731	. ztfi_omp, !! pt
732	. zqfi_omp, !! pq
733	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
734	. zdufi_omp, !! pdu
735	. zdvfi_omp, !! pdv
736	. zdtfi_omp, !! pdt
737	. zdqfi_omp, !! pdq
738	. zdpsrf_omp, !! pdpsrf
739	. tracerdyn) !! tracerdyn <-- utilite ???
740
741	else ! a moduler pour Mars
742	CALL physiq (klon,
743	. llm,
744	. nqtot,
745	. debut_split,
746	. lafin_split,
747	. jD_cur,
748	. jH_cur_split,
749	. zdt_split,
750	. zplev_omp,
751	. zplay_omp,
752	. zpk_omp,
753	. zphi_omp,
754	. zphis_omp,
755	. presnivs_omp,
756	. zufi_omp,
757	. zvfi_omp,
758	. ztfi_omp,
759	. zqfi_omp,
760	. flxwfi_omp,
761	. zdufi_omp,
762	. zdvfi_omp,
763	. zdtfi_omp,
764	. zdqfi_omp,
765	. zdpsrf_omp)
766	endif ! planet_type
767	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
768	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
769	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
770	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
771
772	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
773	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
774	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
775	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
776
777	enddo
778
779	#endif
780	! of #ifdef CPP_PHYS
781
782	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
783	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
784	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
785	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
786
787	c$OMP BARRIER
788
789	do l=1,llm+1
790	do i=1,klon
791	zplev(offset+i,l)=zplev_omp(i,l)
792	enddo
793	enddo
794
795	do l=1,llm
796	do i=1,klon
797	zplay(offset+i,l)=zplay_omp(i,l)
798	enddo
799	enddo
800
801	do l=1,llm
802	do i=1,klon
803	zphi(offset+i,l)=zphi_omp(i,l)
804	enddo
805	enddo
806
807
808	do i=1,klon
809	zphis(offset+i)=zphis_omp(i)
810	enddo
811
812
813	do l=1,llm
814	presnivs(l)=presnivs_omp(l)
815	enddo
816
817	do l=1,llm
818	do i=1,klon
819	zufi(offset+i,l)=zufi_omp(i,l)
820	enddo
821	enddo
822
823	do l=1,llm
824	do i=1,klon
825	zvfi(offset+i,l)=zvfi_omp(i,l)
826	enddo
827	enddo
828
829	do l=1,llm
830	do i=1,klon
831	ztfi(offset+i,l)=ztfi_omp(i,l)
832	enddo
833	enddo
834
835	do iq=1,nqtot
836	do l=1,llm
837	do i=1,klon
838	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
839	enddo
840	enddo
841	enddo
842
843	do l=1,llm
844	do i=1,klon
845	zdufi(offset+i,l)=zdufi_omp(i,l)
846	enddo
847	enddo
848
849	do l=1,llm
850	do i=1,klon
851	zdvfi(offset+i,l)=zdvfi_omp(i,l)
852	enddo
853	enddo
854
855	do l=1,llm
856	do i=1,klon
857	zdtfi(offset+i,l)=zdtfi_omp(i,l)
858	enddo
859	enddo
860
861	do iq=1,nqtot
862	do l=1,llm
863	do i=1,klon
864	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
865	enddo
866	enddo
867	enddo
868
869	do i=1,klon
870	zdpsrf(offset+i)=zdpsrf_omp(i)
871	enddo
872
873
874	klon=klon_mpi
875	500 CONTINUE
876	c$OMP BARRIER
877
878	c$OMP MASTER
879	call stop_timer(timer_physic)
880	c$OMP END MASTER
881
882	IF (using_mpi) THEN
883
884	if (MPI_rank>0) then
885
886	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
887	DO l=1,llm
888	du_send(1:iim,l)=zdufi(1:iim,l)
889	dv_send(1:iim,l)=zdvfi(1:iim,l)
890	ENDDO
891	c$OMP END DO NOWAIT
892
893	c$OMP BARRIER
894	#ifdef CPP_MPI
895	c$OMP MASTER
896	!$OMP CRITICAL (MPI)
897	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
898	& COMM_LMDZ,Req(1),ierr)
899	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
900	& COMM_LMDZ,Req(2),ierr)
901	!$OMP END CRITICAL (MPI)
902	c$OMP END MASTER
903	#endif
904	c$OMP BARRIER
905
906	endif
907
908	if (MPI_rank<MPI_Size-1) then
909	c$OMP BARRIER
910	#ifdef CPP_MPI
911	c$OMP MASTER
912	!$OMP CRITICAL (MPI)
913	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
914	& COMM_LMDZ,Req(3),ierr)
915	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
916	& COMM_LMDZ,Req(4),ierr)
917	!$OMP END CRITICAL (MPI)
918	c$OMP END MASTER
919	#endif
920	endif
921
922	c$OMP BARRIER
923
924
925	#ifdef CPP_MPI
926	c$OMP MASTER
927	!$OMP CRITICAL (MPI)
928	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
929	call MPI_WAITALL(4,Req(1),Status,ierr)
930	else if (MPI_rank>0) then
931	call MPI_WAITALL(2,Req(1),Status,ierr)
932	else if (MPI_rank <MPI_Size-1) then
933	call MPI_WAITALL(2,Req(3),Status,ierr)
934	endif
935	!$OMP END CRITICAL (MPI)
936	c$OMP END MASTER
937	#endif
938
939	c$OMP BARRIER
940
941	ENDIF ! using_mpi
942
943
944	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
945	DO l=1,llm
946
947	zdufi2(1:klon,l)=zdufi(1:klon,l)
948	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
949
950	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
951	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
952
953	pdhfi(:,jj_begin,l)=0
954	pdqfi(:,jj_begin,l,:)=0
955	pdufi(:,jj_begin,l)=0
956	pdvfi(:,jj_begin,l)=0
957
958	if (.not. is_south_pole) then
959	pdhfi(:,jj_end,l)=0
960	pdqfi(:,jj_end,l,:)=0
961	pdufi(:,jj_end,l)=0
962	pdvfi(:,jj_end,l)=0
963	endif
964
965	ENDDO
966	c$OMP END DO NOWAIT
967
968	c$OMP MASTER
969	pdpsfi(:,jj_begin)=0
970	if (.not. is_south_pole) then
971	pdpsfi(:,jj_end)=0
972	endif
973	c$OMP END MASTER
974	c-----------------------------------------------------------------------
975	c transformation des tendances physiques en tendances dynamiques:
976	c ---------------------------------------------------------------
977
978	c tendance sur la pression :
979	c -----------------------------------
980	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
981	c
982	c 62. enthalpie potentielle
983	c ---------------------
984
985	kstart=1
986	kend=klon
987
988	if (is_north_pole) kstart=2
989	if (is_south_pole) kend=klon-1
990
991	! ADAPTATION GCM POUR CP(T)
992	call t2tpot_p(klon,llm,ztfi,zteta,zpk)
993
994
995	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
996	DO l=1,llm
997	!CDIR ON_ADB(index_i)
998	!CDIR ON_ADB(index_j)
999	!cdir NODEP
1000	do ig0=kstart,kend
1001	i=index_i(ig0)
1002	j=index_j(ig0)
1003	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1004	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1005	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1006	if (i==1) then
1007	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
1008	endif
1009	enddo
1010
1011	if (is_north_pole) then
1012	DO i=1,iip1
1013	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
1014	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
1015	enddo
1016	endif
1017
1018	if (is_south_pole) then
1019	DO i=1,iip1
1020	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1021	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
1022	ENDDO
1023	endif
1024	ENDDO
1025	c$OMP END DO NOWAIT
1026
1027	c 62. humidite specifique
1028	c ---------------------
1029	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1030	! DO iq=1,nqtot
1031	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1032	! DO l=1,llm
1033	!!!cdir NODEP
1034	! do ig0=kstart,kend
1035	! i=index_i(ig0)
1036	! j=index_j(ig0)
1037	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1038	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1039	! enddo
1040	!
1041	! if (is_north_pole) then
1042	! do i=1,iip1
1043	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1044	! enddo
1045	! endif
1046	!
1047	! if (is_south_pole) then
1048	! do i=1,iip1
1049	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1050	! enddo
1051	! endif
1052	! ENDDO
1053	!c$OMP END DO NOWAIT
1054	! ENDDO
1055
1056	c 63. traceurs (tous en intensifs)
1057	c ------------
1058	C initialisation des tendances
1059
1060	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1061	DO l=1,llm
1062	pdqfi(:,:,l,:)=0.
1063	ENDDO
1064	c$OMP END DO NOWAIT
1065
1066	C
1067	!cdir NODEP
1068	DO iq=1,nqtot
1069	iiq=niadv(iq)
1070	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1071	DO l=1,llm
1072	!CDIR ON_ADB(index_i)
1073	!CDIR ON_ADB(index_j)
1074	!cdir NODEP
1075	DO ig0=kstart,kend
1076	i=index_i(ig0)
1077	j=index_j(ig0)
1078	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1079	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1080	ENDDO
1081
1082	IF (is_north_pole) then
1083	DO i=1,iip1
1084	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1085	ENDDO
1086	ENDIF
1087
1088	IF (is_south_pole) then
1089	DO i=1,iip1
1090	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1091	ENDDO
1092	ENDIF
1093
1094	ENDDO
1095	c$OMP END DO NOWAIT
1096	ENDDO
1097
1098	c 65. champ u:
1099	c ------------
1100	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1101	DO l=1,llm
1102	!CDIR ON_ADB(index_i)
1103	!CDIR ON_ADB(index_j)
1104	!cdir NODEP
1105	do ig0=kstart,kend
1106	i=index_i(ig0)
1107	j=index_j(ig0)
1108
1109	if (i/=iim) then
1110	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1111	endif
1112
1113	if (i==1) then
1114	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1115	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1116	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1117	endif
1118
1119	enddo
1120
1121	if (is_north_pole) then
1122	DO i=1,iip1
1123	pdufi(i,1,l) = 0.
1124	ENDDO
1125	endif
1126
1127	if (is_south_pole) then
1128	DO i=1,iip1
1129	pdufi(i,jjp1,l) = 0.
1130	ENDDO
1131	endif
1132
1133	ENDDO
1134	c$OMP END DO NOWAIT
1135
1136	c 67. champ v:
1137	c ------------
1138
1139	kstart=1
1140	kend=klon
1141
1142	if (is_north_pole) kstart=2
1143	if (is_south_pole) kend=klon-1-iim
1144
1145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1146	DO l=1,llm
1147	!CDIR ON_ADB(index_i)
1148	!CDIR ON_ADB(index_j)
1149	!cdir NODEP
1150	do ig0=kstart,kend
1151	i=index_i(ig0)
1152	j=index_j(ig0)
1153	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1154	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1155	$ zdvfi2(ig0+iim,l))
1156	$ *cv(i,j)
1157	enddo
1158
1159	ENDDO
1160	c$OMP END DO NOWAIT
1161
1162
1163	c 68. champ v pres des poles:
1164	c ---------------------------
1165	c v = U * cos(long) + V * SIN(long)
1166
1167	if (is_north_pole) then
1168
1169	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1170	DO l=1,llm
1171
1172	DO i=1,iim
1173	pdvfi(i,1,l)=
1174	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1175
1176	pdvfi(i,1,l)=
1177	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1178	ENDDO
1179
1180	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1181
1182	ENDDO
1183	c$OMP END DO NOWAIT
1184
1185	endif
1186
1187	if (is_south_pole) then
1188
1189	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1190	DO l=1,llm
1191
1192	DO i=1,iim
1193	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1194	$ +zdvfi(klon,l)*SIN(rlonv(i))
1195
1196	pdvfi(i,jjm,l)=
1197	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1198	ENDDO
1199
1200	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1201
1202	ENDDO
1203	c$OMP END DO NOWAIT
1204
1205	endif
1206	c-----------------------------------------------------------------------
1207
1208	700 CONTINUE
1209
1210	firstcal = .FALSE.
1211
1212	#else
1213	write(lunout,*)
1214	& "calfis_p: for now can only work with parallel physics"
1215	stop
1216	#endif
1217	! of #ifdef CPP_PHYS
1218	RETURN
1219	END

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

Download in other formats:

Original Format