Context Navigation

calfis_p.F @ 1198

Last change on this file since 1198 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.1 KB

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

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

Download in other formats:

Original Format