Context Navigation

calfis_p.F @ 1356

Last change on this file since 1356 was 1325, checked in by idelkadi, 15 years ago
Decoupage du pas de temps physique en N sous-pas (splitting)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.8 KB

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

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

Download in other formats:

Original Format