Context Navigation

calfis_p.F @ 985

Last change on this file since 985 was 985, checked in by Laurent Fairhead, 16 years ago
Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.9 KB

Line
1	!
2	! $Header$
3	!
4	C
5	C
6	SUBROUTINE calfis_p(nq,
7	$ lafin,
8	$ rdayvrai,
9	$ heure,
10	$ pucov,
11	$ pvcov,
12	$ pteta,
13	$ pq,
14	$ pmasse,
15	$ pps,
16	$ pp,
17	$ ppk,
18	$ pphis,
19	$ pphi,
20	$ pducov,
21	$ pdvcov,
22	$ pdteta,
23	$ pdq,
24	$ flxw,
25	$ clesphy0,
26	$ pdufi,
27	$ pdvfi,
28	$ pdhfi,
29	$ pdqfi,
30	$ pdpsfi)
31	c
32	c Auteur : P. Le Van, F. Hourdin
33	c .........
34	USE dimphy
35	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
36	USE parallel, ONLY : omp_chunk
37	USE mod_interface_dyn_phys
38	USE Write_Field
39	Use Write_field_p
40	USE Times
41	USE IOPHY
42	IMPLICIT NONE
43	c=======================================================================
44	c
45	c 1. rearrangement des tableaux et transformation
46	c variables dynamiques > variables physiques
47	c 2. calcul des termes physiques
48	c 3. retransformation des tendances physiques en tendances dynamiques
49	c
50	c remarques:
51	c ----------
52	c
53	c - les vents sont donnes dans la physique par leurs composantes
54	c naturelles.
55	c - la variable thermodynamique de la physique est une variable
56	c intensive : T
57	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
58	c - les deux seules variables dependant de la geometrie necessaires
59	c pour la physique sont la latitude pour le rayonnement et
60	c l'aire de la maille quand on veut integrer une grandeur
61	c horizontalement.
62	c - les points de la physique sont les points scalaires de la
63	c la dynamique; numerotation:
64	c 1 pour le pole nord
65	c (jjm-1)*iim pour l'interieur du domaine
66	c ngridmx pour le pole sud
67	c ---> ngridmx=2+(jjm-1)*iim
68	c
69	c Input :
70	c -------
71	c ecritphy frequence d'ecriture (en jours)de histphy
72	c pucov covariant zonal velocity
73	c pvcov covariant meridional velocity
74	c pteta potential temperature
75	c pps surface pressure
76	c pmasse masse d'air dans chaque maille
77	c pts surface temperature (K)
78	c callrad clef d'appel au rayonnement
79	c
80	c Output :
81	c --------
82	c pdufi tendency for the natural zonal velocity (ms-1)
83	c pdvfi tendency for the natural meridional velocity
84	c pdhfi tendency for the potential temperature
85	c pdtsfi tendency for the surface temperature
86	c
87	c pdtrad radiative tendencies \ both input
88	c pfluxrad radiative fluxes / and output
89	c
90	c=======================================================================
91	c
92	c-----------------------------------------------------------------------
93	c
94	c 0. Declarations :
95	c ------------------
96
97	#include "dimensions.h"
98	#include "paramet.h"
99	#include "temps.h"
100	#include "advtrac.h"
101
102	INTEGER ngridmx,nq
103	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
104
105	#include "comconst.h"
106	#include "comvert.h"
107	#include "comgeom2.h"
108	#include "control.h"
109	include 'mpif.h'
110
111	c Arguments :
112	c -----------
113	LOGICAL lafin
114	REAL heure
115
116	REAL pvcov(iip1,jjm,llm)
117	REAL pucov(iip1,jjp1,llm)
118	REAL pteta(iip1,jjp1,llm)
119	REAL pmasse(iip1,jjp1,llm)
120	REAL pq(iip1,jjp1,llm,nqmx)
121	REAL pphis(iip1,jjp1)
122	REAL pphi(iip1,jjp1,llm)
123	c
124	REAL pdvcov(iip1,jjm,llm)
125	REAL pducov(iip1,jjp1,llm)
126	REAL pdteta(iip1,jjp1,llm)
127	REAL pdq(iip1,jjp1,llm,nqmx)
128	c
129	REAL pps(iip1,jjp1)
130	REAL pp(iip1,jjp1,llmp1)
131	REAL ppk(iip1,jjp1,llm)
132	c
133	REAL pdvfi(iip1,jjm,llm)
134	REAL pdufi(iip1,jjp1,llm)
135	REAL pdhfi(iip1,jjp1,llm)
136	REAL pdqfi(iip1,jjp1,llm,nqmx)
137	REAL pdpsfi(iip1,jjp1)
138
139	INTEGER longcles
140	PARAMETER ( longcles = 20 )
141	REAL clesphy0( longcles )
142
143
144	c Local variables :
145	c -----------------
146
147	INTEGER i,j,l,ig0,ig,iq,iiq
148	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
149	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
150	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
151	c
152	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
153	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
154	c
155	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
156	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
157	c
158	c REAL,ALLOCATABLE,SAVE :: pvervel(:,:)
159	c
160	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
161	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
162	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
163	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
164
165	c
166	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
167	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
168	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
169	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
170	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
171	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
172	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
173	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
174	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
175	c REAL,ALLOCATABLE,SAVE :: pvervel_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	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
184	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
185	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
186	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp)
187
188	LOGICAL,SAVE :: first_omp=.true.
189	c$OMP THREADPRIVATE(first_omp)
190
191	REAL zsin(iim),zcos(iim),z1(iim)
192	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
193	REAL unskap, pksurcp
194	c
195	cIM diagnostique PVteta, Amip2
196	INTEGER ntetaSTD
197	PARAMETER(ntetaSTD=3)
198	REAL rtetaSTD(ntetaSTD)
199	DATA rtetaSTD/350., 380., 405./
200	REAL PVteta(klon,ntetaSTD)
201
202	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
203
204	REAL SSUM
205
206	LOGICAL firstcal, debut
207	DATA firstcal/.true./
208	SAVE firstcal,debut
209	c$OMP THREADPRIVATE(firstcal,debut)
210	REAL rdayvrai
211
212	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
213	INTEGER :: ierr
214	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
215	INTEGER, dimension(4) :: Req
216	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
217	integer :: k,kstart,kend
218	INTEGER :: offset
219	c
220	c-----------------------------------------------------------------------
221	c
222	c 1. Initialisations :
223	c --------------------
224	c
225
226	klon=klon_mpi
227
228	PVteta(:,:)=0.
229
230	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
231	PRINT*,'STOP dans calfis'
232	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
233	PRINT*,' ngridmx jjm iim '
234	PRINT*,ngridmx,jjm,iim
235	STOP
236	ENDIF
237
238	c-----------------------------------------------------------------------
239	c latitude, longitude et aires des mailles pour la physique:
240	c ----------------------------------------------------------
241
242	c
243	IF ( firstcal ) THEN
244	debut = .TRUE.
245	c$OMP MASTER
246	ALLOCATE(zpsrf(klon))
247	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
248	ALLOCATE(zphi(klon,llm),zphis(klon))
249	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
250	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqmx))
251	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
252	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
253	c ALLOCATE(pvervel(klon,llm))
254	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
255	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqmx))
256	ALLOCATE(zdpsrf(klon))
257	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
258	ALLOCATE(flxwfi(klon,llm))
259	c$OMP END MASTER
260	c$OMP BARRIER
261	ELSE
262	debut = .FALSE.
263	ENDIF
264
265	c
266	c
267	c-----------------------------------------------------------------------
268	c 40. transformation des variables dynamiques en variables physiques:
269	c ---------------------------------------------------------------
270
271	c 41. pressions au sol (en Pascals)
272	c ----------------------------------
273
274	c$OMP MASTER
275	call start_timer(timer_physic)
276	c$OMP END MASTER
277
278	c$OMP MASTER
279	do ig0=1,klon
280	i=index_i(ig0)
281	j=index_j(ig0)
282	zpsrf(ig0)=pps(i,j)
283	enddo
284	c$OMP END MASTER
285
286
287	c 42. pression intercouches :
288	c
289	c -----------------------------------------------------------------
290	c .... zplev definis aux (llm +1) interfaces des couches ....
291	c .... zplay definis aux ( llm ) milieux des couches ....
292	c -----------------------------------------------------------------
293
294	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
295	c
296	unskap = 1./ kappa
297	c
298	c print *,omp_rank,'klon--->',klon
299	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
300	DO l = 1, llmp1
301	do ig0=1,klon
302	i=index_i(ig0)
303	j=index_j(ig0)
304	zplev( ig0,l ) = pp(i,j,l)
305	enddo
306	ENDDO
307	c$OMP END DO NOWAIT
308	c
309	c
310
311	c 43. temperature naturelle (en K) et pressions milieux couches .
312	c ---------------------------------------------------------------
313	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
314	DO l=1,llm
315
316	do ig0=1,klon
317	i=index_i(ig0)
318	j=index_j(ig0)
319	pksurcp = ppk(i,j,l) / cpp
320	zplay(ig0,l) = preff * pksurcp ** unskap
321	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
322	c pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
323	enddo
324
325	ENDDO
326	c$OMP END DO NOWAIT
327
328	c 43.bis traceurs
329	c ---------------
330	c
331
332	DO iq=1,nq
333	iiq=niadv(iq)
334	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
335	DO l=1,llm
336	do ig0=1,klon
337	i=index_i(ig0)
338	j=index_j(ig0)
339	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
340	enddo
341	ENDDO
342	c$OMP END DO NOWAIT
343	ENDDO
344
345	c convergence dynamique pour les traceurs "EAU"
346
347	DO iq=1,2
348	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
349	DO l=1,llm
350	do ig0=1,klon
351	i=index_i(ig0)
352	j=index_j(ig0)
353	c pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
354	enddo
355	ENDDO
356	c$OMP END DO NOWAIT
357	ENDDO
358
359
360
361	c Geopotentiel calcule par rapport a la surface locale:
362	c -----------------------------------------------------
363
364	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
365
366	c$OMP MASTER
367	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
368	c$OMP END MASTER
369	c$OMP BARRIER
370
371	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
372	DO l=1,llm
373	DO ig=1,klon
374	zphi(ig,l)=zphi(ig,l)-zphis(ig)
375	ENDDO
376	ENDDO
377	c$OMP END DO NOWAIT
378
379	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
380	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
381	c de masse est calclue dans advtrac_p.F
382	c
383	cc$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
384	c DO l=1,llm
385	c do ig0=1,klon
386	c i=index_i(ig0)
387	c j=index_j(ig0)
388	c pvervel(ig0,l) = pw(i,j,l)g unsaire(i,j)
389	c enddo
390	c if (is_north_pole) pvervel(1,l)=pw(1,1,l)*g /apoln
391	c if (is_south_pole) pvervel(klon,l)=pw(1,jjp1,l)*g/apols
392	c ENDDO
393	cc$OMP END DO NOWAIT
394
395	c
396	c 45. champ u:
397	c ------------
398
399	kstart=1
400	kend=klon
401
402	if (is_north_pole) kstart=2
403	if (is_south_pole) kend=klon-1
404
405	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
406	DO l=1,llm
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	c pcvgu(ig0,l)= 0.5*( pducov(iim,j,l)/cu(iim,j)
414	c $ + pducov(1,j,l)/cu(1,j) )
415	else
416	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
417	$ + pucov(i,j,l)/cu(i,j) )
418	c pcvgu(ig0,l)= 0.5*( pducov(i-1,j,l)/cu(i-1,j)
419	c $ + pducov(i,j,l)/cu(i,j) )
420	endif
421	enddo
422	ENDDO
423	c$OMP END DO NOWAIT
424	c 46.champ v:
425	c -----------
426	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
427	DO l=1,llm
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	c pcvgv(ig0+i,l)= 0.5 * ( pdvcov(i,j-1,l)/cv(i,j-1)
435	c $ + pdvcov(i,j,l)/cv(i,j) )
436	ENDDO
437	ENDDO
438	c$OMP END DO NOWAIT
439
440	c 47. champs de vents aux pole nord
441	c ------------------------------
442	c U = 1 / pi * integrale [ v * cos(long) * d long ]
443	c V = 1 / pi * integrale [ v * sin(long) * d long ]
444
445	if (is_north_pole) then
446	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
447	DO l=1,llm
448
449	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
450	c z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
451	DO i=2,iim
452	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
453	c z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
454	ENDDO
455
456	DO i=1,iim
457	zcos(i) = COS(rlonv(i))*z1(i)
458	c zcosbis(i)= COS(rlonv(i))*z1bis(i)
459	zsin(i) = SIN(rlonv(i))*z1(i)
460	c zsinbis(i)= SIN(rlonv(i))*z1bis(i)
461	ENDDO
462
463	zufi(1,l) = SSUM(iim,zcos,1)/pi
464	c pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
465	zvfi(1,l) = SSUM(iim,zsin,1)/pi
466	c pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
467
468	ENDDO
469	c$OMP END DO NOWAIT
470	endif
471
472
473	c 48. champs de vents aux pole sud:
474	c ---------------------------------
475	c U = 1 / pi * integrale [ v * cos(long) * d long ]
476	c V = 1 / pi * integrale [ v * sin(long) * d long ]
477
478	if (is_south_pole) then
479	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
480	DO l=1,llm
481
482	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
483	c z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
484	DO i=2,iim
485	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
486	c z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
487	ENDDO
488
489	DO i=1,iim
490	zcos(i) = COS(rlonv(i))*z1(i)
491	c zcosbis(i) = COS(rlonv(i))*z1bis(i)
492	zsin(i) = SIN(rlonv(i))*z1(i)
493	c zsinbis(i) = SIN(rlonv(i))*z1bis(i)
494	ENDDO
495
496	zufi(klon,l) = SSUM(iim,zcos,1)/pi
497	c pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
498	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
499	c pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
500
501	ENDDO
502	c$OMP END DO NOWAIT
503	endif
504
505
506	IF (is_sequential) THEN
507	c
508	cIM calcul PV a teta=350, 380, 405K
509	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
510	$ ztfi,zplay,zplev,
511	$ ntetaSTD,rtetaSTD,PVteta)
512	c
513	ENDIF
514
515	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
516	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
517
518	c-----------------------------------------------------------------------
519	c Appel de la physique:
520	c ---------------------
521
522	cc$OMP PARALLEL DEFAULT(NONE)
523	cc$OMP+ PRIVATE(i,l,offset,iq)
524	cc$OMP+ SHARED(klon_omp_nb,nq,klon_omp_begin,
525	cc$OMP+ debut,lafin,rdayvrai,heure,dtphys,zplev,zplay,
526	cc$OMP+ zphi,zphis,presnivs,clesphy0,zufi,zvfi,ztfi,
527	cc$OMP+ zqfi,pvervel,zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
528
529	c PRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
530	c c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
531	c c$OMP+ zqfi_omp,pvervel_omp,zdufi_omp,zdvfi_omp,
532	c c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp)
533
534	c$OMP BARRIER
535	if (first_omp) then
536	klon=klon_omp
537
538	allocate(zplev_omp(klon,llm+1))
539	allocate(zplay_omp(klon,llm))
540	allocate(zphi_omp(klon,llm))
541	allocate(zphis_omp(klon))
542	allocate(presnivs_omp(llm))
543	allocate(zufi_omp(klon,llm))
544	allocate(zvfi_omp(klon,llm))
545	allocate(ztfi_omp(klon,llm))
546	allocate(zqfi_omp(klon,llm,nq))
547	c allocate(pvervel_omp(klon,llm))
548	allocate(zdufi_omp(klon,llm))
549	allocate(zdvfi_omp(klon,llm))
550	allocate(zdtfi_omp(klon,llm))
551	allocate(zdqfi_omp(klon,llm,nq))
552	allocate(zdpsrf_omp(klon))
553	allocate(flxwfi_omp(klon,llm))
554	first_omp=.false.
555	endif
556
557
558	klon=klon_omp
559	offset=klon_omp_begin-1
560
561	do l=1,llm+1
562	do i=1,klon
563	zplev_omp(i,l)=zplev(offset+i,l)
564	enddo
565	enddo
566
567	do l=1,llm
568	do i=1,klon
569	zplay_omp(i,l)=zplay(offset+i,l)
570	enddo
571	enddo
572
573	do l=1,llm
574	do i=1,klon
575	zphi_omp(i,l)=zphi(offset+i,l)
576	enddo
577	enddo
578
579	do i=1,klon
580	zphis_omp(i)=zphis(offset+i)
581	enddo
582
583
584	do l=1,llm
585	presnivs_omp(l)=presnivs(l)
586	enddo
587
588	do l=1,llm
589	do i=1,klon
590	zufi_omp(i,l)=zufi(offset+i,l)
591	enddo
592	enddo
593
594	do l=1,llm
595	do i=1,klon
596	zvfi_omp(i,l)=zvfi(offset+i,l)
597	enddo
598	enddo
599
600	do l=1,llm
601	do i=1,klon
602	ztfi_omp(i,l)=ztfi(offset+i,l)
603	enddo
604	enddo
605
606	do iq=1,nq
607	do l=1,llm
608	do i=1,klon
609	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
610	enddo
611	enddo
612	enddo
613
614	c do l=1,llm
615	c do i=1,klon
616	c pvervel_omp(i,l)=pvervel(offset+i,l)
617	c enddo
618	c enddo
619
620	do l=1,llm
621	do i=1,klon
622	zdufi_omp(i,l)=zdufi(offset+i,l)
623	enddo
624	enddo
625
626	do l=1,llm
627	do i=1,klon
628	zdvfi_omp(i,l)=zdvfi(offset+i,l)
629	enddo
630	enddo
631
632	do l=1,llm
633	do i=1,klon
634	zdtfi_omp(i,l)=zdtfi(offset+i,l)
635	enddo
636	enddo
637
638	do iq=1,nq
639	do l=1,llm
640	do i=1,klon
641	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
642	enddo
643	enddo
644	enddo
645
646	do i=1,klon
647	zdpsrf_omp(i)=zdpsrf(offset+i)
648	enddo
649
650	do l=1,llm
651	do i=1,klon
652	flxwfi_omp(i,l)=flxwfi(offset+i,l)
653	enddo
654	enddo
655
656	c$OMP BARRIER
657	cym call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
658
659	CALL physiq (klon,
660	. llm,
661	. nq,
662	. debut,
663	. lafin,
664	. rdayvrai,
665	. heure,
666	. dtphys,
667	. zplev_omp,
668	. zplay_omp,
669	. zphi_omp,
670	. zphis_omp,
671	. presnivs_omp,
672	. clesphy0,
673	. zufi_omp,
674	. zvfi_omp,
675	. ztfi_omp,
676	. zqfi_omp,
677	c . pvervel_omp,
678	c#ifdef INCA
679	. flxwfi_omp,
680	c#endif
681	. zdufi_omp,
682	. zdvfi_omp,
683	. zdtfi_omp,
684	. zdqfi_omp,
685	. zdpsrf_omp,
686	cIM diagnostique PVteta, Amip2
687	. pducov,
688	. PVteta)
689
690	cym call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
691
692	c$OMP BARRIER
693
694	do l=1,llm+1
695	do i=1,klon
696	zplev(offset+i,l)=zplev_omp(i,l)
697	enddo
698	enddo
699
700	do l=1,llm
701	do i=1,klon
702	zplay(offset+i,l)=zplay_omp(i,l)
703	enddo
704	enddo
705
706	do l=1,llm
707	do i=1,klon
708	zphi(offset+i,l)=zphi_omp(i,l)
709	enddo
710	enddo
711
712
713	do i=1,klon
714	zphis(offset+i)=zphis_omp(i)
715	enddo
716
717
718	do l=1,llm
719	presnivs(l)=presnivs_omp(l)
720	enddo
721
722	do l=1,llm
723	do i=1,klon
724	zufi(offset+i,l)=zufi_omp(i,l)
725	enddo
726	enddo
727
728	do l=1,llm
729	do i=1,klon
730	zvfi(offset+i,l)=zvfi_omp(i,l)
731	enddo
732	enddo
733
734	do l=1,llm
735	do i=1,klon
736	ztfi(offset+i,l)=ztfi_omp(i,l)
737	enddo
738	enddo
739
740	do iq=1,nq
741	do l=1,llm
742	do i=1,klon
743	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
744	enddo
745	enddo
746	enddo
747
748	c do l=1,llm
749	c do i=1,klon
750	c pvervel(offset+i,l)=pvervel_omp(i,l)
751	c enddo
752	c enddo
753
754	do l=1,llm
755	do i=1,klon
756	zdufi(offset+i,l)=zdufi_omp(i,l)
757	enddo
758	enddo
759
760	do l=1,llm
761	do i=1,klon
762	zdvfi(offset+i,l)=zdvfi_omp(i,l)
763	enddo
764	enddo
765
766	do l=1,llm
767	do i=1,klon
768	zdtfi(offset+i,l)=zdtfi_omp(i,l)
769	enddo
770	enddo
771
772	do iq=1,nq
773	do l=1,llm
774	do i=1,klon
775	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
776	enddo
777	enddo
778	enddo
779
780	do i=1,klon
781	zdpsrf(offset+i)=zdpsrf_omp(i)
782	enddo
783
784
785	cc$OMP END PARALLEL
786	klon=klon_mpi
787	500 CONTINUE
788	c$OMP BARRIER
789
790	c$OMP MASTER
791	cym call WriteField_phy('zdtfi',zdtfi(:,:),llm)
792	call stop_timer(timer_physic)
793	c$OMP END MASTER
794
795	if (MPI_rank>0) then
796
797	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
798	DO l=1,llm
799	du_send(1:iim,l)=zdufi(1:iim,l)
800	dv_send(1:iim,l)=zdvfi(1:iim,l)
801	ENDDO
802	c$OMP END DO NOWAIT
803
804	c$OMP BARRIER
805	c$OMP MASTER
806	!$OMP CRITICAL (MPI)
807	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
808	& COMM_LMDZ,Req(1),ierr)
809	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
810	& COMM_LMDZ,Req(2),ierr)
811	!$OMP END CRITICAL (MPI)
812	c$OMP END MASTER
813	c$OMP BARRIER
814
815	endif
816
817	if (MPI_rank<MPI_Size-1) then
818	c$OMP BARRIER
819	c$OMP MASTER
820	!$OMP CRITICAL (MPI)
821	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
822	& COMM_LMDZ,Req(3),ierr)
823	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
824	& COMM_LMDZ,Req(4),ierr)
825	!$OMP END CRITICAL (MPI)
826	c$OMP END MASTER
827	c$OMP BARRIER
828	endif
829
830	c$OMP BARRIER
831	c$OMP MASTER
832	!$OMP CRITICAL (MPI)
833	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
834	call MPI_WAITALL(4,Req(1),Status,ierr)
835	else if (MPI_rank>0) then
836	call MPI_WAITALL(2,Req(1),Status,ierr)
837	else if (MPI_rank <MPI_Size-1) then
838	call MPI_WAITALL(2,Req(3),Status,ierr)
839	endif
840	!$OMP END CRITICAL (MPI)
841	c$OMP END MASTER
842	c$OMP BARRIER
843
844	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
845	DO l=1,llm
846
847	zdufi2(1:klon,l)=zdufi(1:klon,l)
848	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
849
850	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
851	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
852
853	pdhfi(:,jj_begin,l)=0
854	pdqfi(:,jj_begin,l,:)=0
855	pdufi(:,jj_begin,l)=0
856	pdvfi(:,jj_begin,l)=0
857
858	if (.not. is_south_pole) then
859	pdhfi(:,jj_end,l)=0
860	pdqfi(:,jj_end,l,:)=0
861	pdufi(:,jj_end,l)=0
862	pdvfi(:,jj_end,l)=0
863	endif
864
865	ENDDO
866	c$OMP END DO NOWAIT
867
868	c$OMP MASTER
869	pdpsfi(:,jj_begin)=0
870	if (.not. is_south_pole) then
871	pdpsfi(:,jj_end)=0
872	endif
873	c$OMP END MASTER
874	c-----------------------------------------------------------------------
875	c transformation des tendances physiques en tendances dynamiques:
876	c ---------------------------------------------------------------
877
878	c tendance sur la pression :
879	c -----------------------------------
880	c$OMP MASTER
881	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
882	c$OMP END MASTER
883	c
884	c 62. enthalpie potentielle
885	c ---------------------
886
887	kstart=1
888	kend=klon
889
890	if (is_north_pole) kstart=2
891	if (is_south_pole) kend=klon-1
892
893	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
894	DO l=1,llm
895
896	!!cdir NODEP
897	do ig0=kstart,kend
898	i=index_i(ig0)
899	j=index_j(ig0)
900	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
901	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
902	enddo
903
904	if (is_north_pole) then
905	DO i=1,iip1
906	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
907	enddo
908	endif
909
910	if (is_south_pole) then
911	DO i=1,iip1
912	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
913	ENDDO
914	endif
915	ENDDO
916	c$OMP END DO NOWAIT
917
918	c 62. humidite specifique
919	c ---------------------
920
921	DO iq=1,nqmx
922	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
923	DO l=1,llm
924	!!cdir NODEP
925	do ig0=kstart,kend
926	i=index_i(ig0)
927	j=index_j(ig0)
928	pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
929	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
930	enddo
931
932	if (is_north_pole) then
933	do i=1,iip1
934	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
935	enddo
936	endif
937
938	if (is_south_pole) then
939	do i=1,iip1
940	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
941	enddo
942	endif
943
944	ENDDO
945	c$OMP END DO NOWAIT
946	ENDDO
947
948	c 63. traceurs
949	c ------------
950	C initialisation des tendances
951
952	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
953	DO l=1,llm
954	pdqfi(:,:,l,:)=0.
955	ENDDO
956	c$OMP END DO NOWAIT
957
958	C
959
960	DO iq=1,nq
961	iiq=niadv(iq)
962	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
963	DO l=1,llm
964
965	!!cdir NODEP
966	DO ig0=kstart,kend
967	i=index_i(ig0)
968	j=index_j(ig0)
969	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
970	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
971	ENDDO
972
973	IF (is_north_pole) then
974	DO i=1,iip1
975	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
976	ENDDO
977	ENDIF
978
979	IF (is_south_pole) then
980	DO i=1,iip1
981	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
982	ENDDO
983	ENDIF
984
985	ENDDO
986	c$OMP END DO NOWAIT
987	ENDDO
988
989	c 65. champ u:
990	c ------------
991	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
992	DO l=1,llm
993	!!cdir NODEP
994	do ig0=kstart,kend
995	i=index_i(ig0)
996	j=index_j(ig0)
997
998	if (i/=iim) then
999	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1000	endif
1001
1002	if (i==1) then
1003	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1004	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1005	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1006	endif
1007
1008	enddo
1009
1010	if (is_north_pole) then
1011	DO i=1,iip1
1012	pdufi(i,1,l) = 0.
1013	ENDDO
1014	endif
1015
1016	if (is_south_pole) then
1017	DO i=1,iip1
1018	pdufi(i,jjp1,l) = 0.
1019	ENDDO
1020	endif
1021
1022	ENDDO
1023	c$OMP END DO NOWAIT
1024
1025	c 67. champ v:
1026	c ------------
1027
1028	kstart=1
1029	kend=klon
1030
1031	if (is_north_pole) kstart=2
1032	if (is_south_pole) kend=klon-1-iim
1033
1034	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1035	DO l=1,llm
1036	!!cdir NODEP
1037	do ig0=kstart,kend
1038	i=index_i(ig0)
1039	j=index_j(ig0)
1040	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1041	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1042	$ zdvfi2(ig0+iim,l))
1043	$ *cv(i,j)
1044	enddo
1045
1046	ENDDO
1047	c$OMP END DO NOWAIT
1048
1049
1050	c 68. champ v pres des poles:
1051	c ---------------------------
1052	c v = U * cos(long) + V * SIN(long)
1053
1054	if (is_north_pole) then
1055
1056	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1057	DO l=1,llm
1058
1059	DO i=1,iim
1060	pdvfi(i,1,l)=
1061	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1062
1063	pdvfi(i,1,l)=
1064	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1065	ENDDO
1066
1067	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1068
1069	ENDDO
1070	c$OMP END DO NOWAIT
1071
1072	endif
1073
1074	if (is_south_pole) then
1075
1076	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1077	DO l=1,llm
1078
1079	DO i=1,iim
1080	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1081	$ +zdvfi(klon,l)*SIN(rlonv(i))
1082
1083	pdvfi(i,jjm,l)=
1084	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1085	ENDDO
1086
1087	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1088
1089	ENDDO
1090	c$OMP END DO NOWAIT
1091
1092	endif
1093	c-----------------------------------------------------------------------
1094
1095	700 CONTINUE
1096
1097	firstcal = .FALSE.
1098
1099	RETURN
1100	END

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

Download in other formats:

Original Format