Context Navigation

calfis_p.F @ 1208

Last change on this file since 1208 was 1207, checked in by Laurent Fairhead, 15 years ago
Probleme sur la forme des commentaires laisses par emacs, openmp pourrait mal les comprendre L'argument heure disparait de la liste des arguments d'appel a calfis
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.0 KB

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

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

Download in other formats:

Original Format