Context Navigation

calfis_p.F @ 8

Last change on this file since 8 was 8, checked in by emillour, 14 years ago

Debut de mise a jour de la dynamique parallele par rapport aux modifs dans la partie sequentielle.

Mais NON TESTE , car pas (encore) possibilite de compiler et faire tourner cas simple (type newtonien sans physique).

Voir commit_v8.log pour les details.

Ehouarn

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

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

Download in other formats:

Original Format