Context Navigation

calfis_p.F @ 4057

Last change on this file since 4057 was 1407, checked in by Laurent Fairhead, 14 years ago

Some extraneous prints deleted for optimisation purposes

Des prints supprimés pour optimisation

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 1407 2010-07-07 10:31:52Z dcugnet $
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 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
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	REAL pdq(iip1,jjp1,llm,nqtot)
132	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
133	c
134	REAL pps(iip1,jjp1)
135	REAL pp(iip1,jjp1,llmp1)
136	REAL ppk(iip1,jjp1,llm)
137	c
138	REAL pdvfi(iip1,jjm,llm)
139	REAL pdufi(iip1,jjp1,llm)
140	REAL pdhfi(iip1,jjp1,llm)
141	REAL pdqfi(iip1,jjp1,llm,nqtot)
142	REAL pdpsfi(iip1,jjp1)
143
144	INTEGER longcles
145	PARAMETER ( longcles = 20 )
146	REAL clesphy0( longcles )
147
148	#ifdef CPP_EARTH
149	c Local variables :
150	c -----------------
151
152	INTEGER i,j,l,ig0,ig,iq,iiq
153	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
154	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
155	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
156	c
157	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
158	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
159	c
160	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
161	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
162	c
163	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
164	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
165	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
166	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
167
168	c
169	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
170	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
171	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
172	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
173	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
174	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
175	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
176	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
177	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
178	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
180	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
181	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
182	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
183	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
184
185	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186	! Introduction du splitting (FH)
187	! Question pour Yann :
188	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
189	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
190	! soit allocatable (plutot par exemple que de passer une dimension
191	! dépendant du process en argument des routines) et que, du coup,
192	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
193	! Tu confirmes ?
194	! J'ai suivi le même principe pour les zdufic_omp
195	! Mais c'est surement bien que tu controles.
196	!
197
198	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
199	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
200	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
201	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
202	REAL jH_cur_split,zdt_split
203	LOGICAL debut_split,lafin_split
204	INTEGER isplit
205	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
206
207	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
208	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
209	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
210	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
211	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
212
213	LOGICAL,SAVE :: first_omp=.true.
214	c$OMP THREADPRIVATE(first_omp)
215
216	REAL zsin(iim),zcos(iim),z1(iim)
217	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
218	REAL unskap, pksurcp
219	c
220	cIM diagnostique PVteta, Amip2
221	INTEGER ntetaSTD
222	PARAMETER(ntetaSTD=3)
223	REAL rtetaSTD(ntetaSTD)
224	DATA rtetaSTD/350., 380., 405./
225	REAL PVteta(klon,ntetaSTD)
226
227
228	REAL SSUM
229
230	LOGICAL firstcal, debut
231	DATA firstcal/.true./
232	SAVE firstcal,debut
233	c$OMP THREADPRIVATE(firstcal,debut)
234
235	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
236	INTEGER :: ierr
237	#ifdef CPP_MPI
238	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
239	#else
240	INTEGER,dimension(1,4) :: Status
241	#endif
242	INTEGER, dimension(4) :: Req
243	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
244	integer :: k,kstart,kend
245	INTEGER :: offset
246	c
247	c-----------------------------------------------------------------------
248	c
249	c 1. Initialisations :
250	c --------------------
251	c
252
253	klon=klon_mpi
254
255	PVteta(:,:)=0.
256
257	c
258	IF ( firstcal ) THEN
259	debut = .TRUE.
260	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
261	write(lunout,*) 'STOP dans calfis'
262	write(lunout,*)
263	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
264	write(lunout,*) ' ngridmx jjm iim '
265	write(lunout,*) ngridmx,jjm,iim
266	STOP
267	ENDIF
268	c$OMP MASTER
269	ALLOCATE(zpsrf(klon))
270	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
271	ALLOCATE(zphi(klon,llm),zphis(klon))
272	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
273	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
274	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
275	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
276	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
277	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
278	ALLOCATE(zdpsrf(klon))
279	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
280	ALLOCATE(flxwfi(klon,llm))
281	c$OMP END MASTER
282	c$OMP BARRIER
283	ELSE
284	debut = .FALSE.
285	ENDIF
286
287	c
288	c
289	c-----------------------------------------------------------------------
290	c 40. transformation des variables dynamiques en variables physiques:
291	c ---------------------------------------------------------------
292
293	c 41. pressions au sol (en Pascals)
294	c ----------------------------------
295
296	c$OMP MASTER
297	call start_timer(timer_physic)
298	c$OMP END MASTER
299
300	c$OMP MASTER
301	!CDIR ON_ADB(index_i)
302	!CDIR ON_ADB(index_j)
303	do ig0=1,klon
304	i=index_i(ig0)
305	j=index_j(ig0)
306	zpsrf(ig0)=pps(i,j)
307	enddo
308	c$OMP END MASTER
309
310
311	c 42. pression intercouches :
312	c
313	c -----------------------------------------------------------------
314	c .... zplev definis aux (llm +1) interfaces des couches ....
315	c .... zplay definis aux ( llm ) milieux des couches ....
316	c -----------------------------------------------------------------
317
318	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
319	c
320	unskap = 1./ kappa
321	c
322	c print *,omp_rank,'klon--->',klon
323	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
324	DO l = 1, llmp1
325	!CDIR ON_ADB(index_i)
326	!CDIR ON_ADB(index_j)
327	do ig0=1,klon
328	i=index_i(ig0)
329	j=index_j(ig0)
330	zplev( ig0,l ) = pp(i,j,l)
331	enddo
332	ENDDO
333	c$OMP END DO NOWAIT
334	c
335	c
336
337	c 43. temperature naturelle (en K) et pressions milieux couches .
338	c ---------------------------------------------------------------
339	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
340	DO l=1,llm
341	!CDIR ON_ADB(index_i)
342	!CDIR ON_ADB(index_j)
343	do ig0=1,klon
344	i=index_i(ig0)
345	j=index_j(ig0)
346	pksurcp = ppk(i,j,l) / cpp
347	zplay(ig0,l) = preff * pksurcp ** unskap
348	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
349	enddo
350
351	ENDDO
352	c$OMP END DO NOWAIT
353
354	c 43.bis traceurs
355	c ---------------
356	c
357
358	DO iq=1,nqtot
359	iiq=niadv(iq)
360	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
361	DO l=1,llm
362	!CDIR ON_ADB(index_i)
363	!CDIR ON_ADB(index_j)
364	do ig0=1,klon
365	i=index_i(ig0)
366	j=index_j(ig0)
367	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
368	enddo
369	ENDDO
370	c$OMP END DO NOWAIT
371	ENDDO
372
373
374	c Geopotentiel calcule par rapport a la surface locale:
375	c -----------------------------------------------------
376
377	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
378
379	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
380
381	c$OMP BARRIER
382
383	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
384	DO l=1,llm
385	DO ig=1,klon
386	zphi(ig,l)=zphi(ig,l)-zphis(ig)
387	ENDDO
388	ENDDO
389	c$OMP END DO NOWAIT
390
391
392	c
393	c 45. champ u:
394	c ------------
395
396	kstart=1
397	kend=klon
398
399	if (is_north_pole) kstart=2
400	if (is_south_pole) kend=klon-1
401
402	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
403	DO l=1,llm
404	!CDIR ON_ADB(index_i)
405	!CDIR ON_ADB(index_j)
406	!CDIR SPARSE
407	do ig0=kstart,kend
408	i=index_i(ig0)
409	j=index_j(ig0)
410	if (i==1) then
411	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
412	$ + pucov(1,j,l)/cu(1,j) )
413	else
414	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
415	$ + pucov(i,j,l)/cu(i,j) )
416	endif
417	enddo
418	ENDDO
419	c$OMP END DO NOWAIT
420
421	c 46.champ v:
422	c -----------
423
424	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
425	DO l=1,llm
426	!CDIR ON_ADB(index_i)
427	!CDIR ON_ADB(index_j)
428	DO ig0=kstart,kend
429	i=index_i(ig0)
430	j=index_j(ig0)
431	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
432	$ + pvcov(i,j,l)/cv(i,j) )
433
434	ENDDO
435	ENDDO
436	c$OMP END DO NOWAIT
437
438	c 47. champs de vents aux pole nord
439	c ------------------------------
440	c U = 1 / pi * integrale [ v * cos(long) * d long ]
441	c V = 1 / pi * integrale [ v * sin(long) * d long ]
442
443	if (is_north_pole) then
444	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
445	DO l=1,llm
446
447	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
448	DO i=2,iim
449	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
450	ENDDO
451
452	DO i=1,iim
453	zcos(i) = COS(rlonv(i))*z1(i)
454	zsin(i) = SIN(rlonv(i))*z1(i)
455	ENDDO
456
457	zufi(1,l) = SSUM(iim,zcos,1)/pi
458	zvfi(1,l) = SSUM(iim,zsin,1)/pi
459
460	ENDDO
461	c$OMP END DO NOWAIT
462	endif
463
464
465	c 48. champs de vents aux pole sud:
466	c ---------------------------------
467	c U = 1 / pi * integrale [ v * cos(long) * d long ]
468	c V = 1 / pi * integrale [ v * sin(long) * d long ]
469
470	if (is_south_pole) then
471	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
472	DO l=1,llm
473
474	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
475	DO i=2,iim
476	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
477	ENDDO
478
479	DO i=1,iim
480	zcos(i) = COS(rlonv(i))*z1(i)
481	zsin(i) = SIN(rlonv(i))*z1(i)
482	ENDDO
483
484	zufi(klon,l) = SSUM(iim,zcos,1)/pi
485	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
486	ENDDO
487	c$OMP END DO NOWAIT
488	endif
489
490
491	IF (is_sequential) THEN
492	c
493	cIM calcul PV a teta=350, 380, 405K
494	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
495	$ ztfi,zplay,zplev,
496	$ ntetaSTD,rtetaSTD,PVteta)
497	c
498	ENDIF
499
500	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
501	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
502
503	c-----------------------------------------------------------------------
504	c Appel de la physique:
505	c ---------------------
506
507
508	c$OMP BARRIER
509	if (first_omp) then
510	klon=klon_omp
511
512	allocate(zplev_omp(klon,llm+1))
513	allocate(zplay_omp(klon,llm))
514	allocate(zphi_omp(klon,llm))
515	allocate(zphis_omp(klon))
516	allocate(presnivs_omp(llm))
517	allocate(zufi_omp(klon,llm))
518	allocate(zvfi_omp(klon,llm))
519	allocate(ztfi_omp(klon,llm))
520	allocate(zqfi_omp(klon,llm,nqtot))
521	allocate(zdufi_omp(klon,llm))
522	allocate(zdvfi_omp(klon,llm))
523	allocate(zdtfi_omp(klon,llm))
524	allocate(zdqfi_omp(klon,llm,nqtot))
525	allocate(zdufic_omp(klon,llm))
526	allocate(zdvfic_omp(klon,llm))
527	allocate(zdtfic_omp(klon,llm))
528	allocate(zdqfic_omp(klon,llm,nqtot))
529	allocate(zdpsrf_omp(klon))
530	allocate(flxwfi_omp(klon,llm))
531	first_omp=.false.
532	endif
533
534
535	klon=klon_omp
536	offset=klon_omp_begin-1
537
538	do l=1,llm+1
539	do i=1,klon
540	zplev_omp(i,l)=zplev(offset+i,l)
541	enddo
542	enddo
543
544	do l=1,llm
545	do i=1,klon
546	zplay_omp(i,l)=zplay(offset+i,l)
547	enddo
548	enddo
549
550	do l=1,llm
551	do i=1,klon
552	zphi_omp(i,l)=zphi(offset+i,l)
553	enddo
554	enddo
555
556	do i=1,klon
557	zphis_omp(i)=zphis(offset+i)
558	enddo
559
560
561	do l=1,llm
562	presnivs_omp(l)=presnivs(l)
563	enddo
564
565	do l=1,llm
566	do i=1,klon
567	zufi_omp(i,l)=zufi(offset+i,l)
568	enddo
569	enddo
570
571	do l=1,llm
572	do i=1,klon
573	zvfi_omp(i,l)=zvfi(offset+i,l)
574	enddo
575	enddo
576
577	do l=1,llm
578	do i=1,klon
579	ztfi_omp(i,l)=ztfi(offset+i,l)
580	enddo
581	enddo
582
583	do iq=1,nqtot
584	do l=1,llm
585	do i=1,klon
586	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
587	enddo
588	enddo
589	enddo
590
591	do l=1,llm
592	do i=1,klon
593	zdufi_omp(i,l)=zdufi(offset+i,l)
594	enddo
595	enddo
596
597	do l=1,llm
598	do i=1,klon
599	zdvfi_omp(i,l)=zdvfi(offset+i,l)
600	enddo
601	enddo
602
603	do l=1,llm
604	do i=1,klon
605	zdtfi_omp(i,l)=zdtfi(offset+i,l)
606	enddo
607	enddo
608
609	do iq=1,nqtot
610	do l=1,llm
611	do i=1,klon
612	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
613	enddo
614	enddo
615	enddo
616
617	do i=1,klon
618	zdpsrf_omp(i)=zdpsrf(offset+i)
619	enddo
620
621	do l=1,llm
622	do i=1,klon
623	flxwfi_omp(i,l)=flxwfi(offset+i,l)
624	enddo
625	enddo
626
627	c$OMP BARRIER
628
629	if (planet_type=="earth") then
630	#ifdef CPP_EARTH
631
632	!$OMP MASTER
633	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
634	!$OMP END MASTER
635	zdt_split=dtphys/nsplit_phys
636	zdufic_omp(:,:)=0.
637	zdvfic_omp(:,:)=0.
638	zdtfic_omp(:,:)=0.
639	zdqfic_omp(:,:,:)=0.
640
641	do isplit=1,nsplit_phys
642
643	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
644	debut_split=debut.and.isplit==1
645	lafin_split=lafin.and.isplit==nsplit_phys
646
647
648	CALL physiq (klon,
649	. llm,
650	. debut_split,
651	. lafin_split,
652	. jD_cur,
653	. jH_cur_split,
654	. zdt_split,
655	. zplev_omp,
656	. zplay_omp,
657	. zphi_omp,
658	. zphis_omp,
659	. presnivs_omp,
660	. clesphy0,
661	. zufi_omp,
662	. zvfi_omp,
663	. ztfi_omp,
664	. zqfi_omp,
665	c#ifdef INCA
666	. flxwfi_omp,
667	c#endif
668	. zdufi_omp,
669	. zdvfi_omp,
670	. zdtfi_omp,
671	. zdqfi_omp,
672	. zdpsrf_omp,
673	cIM diagnostique PVteta, Amip2
674	. pducov,
675	. PVteta)
676
677	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
678	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
679	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
680	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
681
682	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
683	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
684	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
685	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
686
687	enddo
688
689	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
690	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
691	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
692	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
693
694	#endif
695	endif !of if (planet_type=="earth")
696	c$OMP BARRIER
697
698	do l=1,llm+1
699	do i=1,klon
700	zplev(offset+i,l)=zplev_omp(i,l)
701	enddo
702	enddo
703
704	do l=1,llm
705	do i=1,klon
706	zplay(offset+i,l)=zplay_omp(i,l)
707	enddo
708	enddo
709
710	do l=1,llm
711	do i=1,klon
712	zphi(offset+i,l)=zphi_omp(i,l)
713	enddo
714	enddo
715
716
717	do i=1,klon
718	zphis(offset+i)=zphis_omp(i)
719	enddo
720
721
722	do l=1,llm
723	presnivs(l)=presnivs_omp(l)
724	enddo
725
726	do l=1,llm
727	do i=1,klon
728	zufi(offset+i,l)=zufi_omp(i,l)
729	enddo
730	enddo
731
732	do l=1,llm
733	do i=1,klon
734	zvfi(offset+i,l)=zvfi_omp(i,l)
735	enddo
736	enddo
737
738	do l=1,llm
739	do i=1,klon
740	ztfi(offset+i,l)=ztfi_omp(i,l)
741	enddo
742	enddo
743
744	do iq=1,nqtot
745	do l=1,llm
746	do i=1,klon
747	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
748	enddo
749	enddo
750	enddo
751
752	do l=1,llm
753	do i=1,klon
754	zdufi(offset+i,l)=zdufi_omp(i,l)
755	enddo
756	enddo
757
758	do l=1,llm
759	do i=1,klon
760	zdvfi(offset+i,l)=zdvfi_omp(i,l)
761	enddo
762	enddo
763
764	do l=1,llm
765	do i=1,klon
766	zdtfi(offset+i,l)=zdtfi_omp(i,l)
767	enddo
768	enddo
769
770	do iq=1,nqtot
771	do l=1,llm
772	do i=1,klon
773	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
774	enddo
775	enddo
776	enddo
777
778	do i=1,klon
779	zdpsrf(offset+i)=zdpsrf_omp(i)
780	enddo
781
782
783	klon=klon_mpi
784	500 CONTINUE
785	c$OMP BARRIER
786
787	c$OMP MASTER
788	call stop_timer(timer_physic)
789	c$OMP END MASTER
790
791	IF (using_mpi) THEN
792
793	if (MPI_rank>0) then
794
795	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
796	DO l=1,llm
797	du_send(1:iim,l)=zdufi(1:iim,l)
798	dv_send(1:iim,l)=zdvfi(1:iim,l)
799	ENDDO
800	c$OMP END DO NOWAIT
801
802	c$OMP BARRIER
803	#ifdef CPP_MPI
804	c$OMP MASTER
805	!$OMP CRITICAL (MPI)
806	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
807	& COMM_LMDZ,Req(1),ierr)
808	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
809	& COMM_LMDZ,Req(2),ierr)
810	!$OMP END CRITICAL (MPI)
811	c$OMP END MASTER
812	#endif
813	c$OMP BARRIER
814
815	endif
816
817	if (MPI_rank<MPI_Size-1) then
818	c$OMP BARRIER
819	#ifdef CPP_MPI
820	c$OMP MASTER
821	!$OMP CRITICAL (MPI)
822	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
823	& COMM_LMDZ,Req(3),ierr)
824	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
825	& COMM_LMDZ,Req(4),ierr)
826	!$OMP END CRITICAL (MPI)
827	c$OMP END MASTER
828	#endif
829	endif
830
831	c$OMP BARRIER
832
833
834	#ifdef CPP_MPI
835	c$OMP MASTER
836	!$OMP CRITICAL (MPI)
837	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
838	call MPI_WAITALL(4,Req(1),Status,ierr)
839	else if (MPI_rank>0) then
840	call MPI_WAITALL(2,Req(1),Status,ierr)
841	else if (MPI_rank <MPI_Size-1) then
842	call MPI_WAITALL(2,Req(3),Status,ierr)
843	endif
844	!$OMP END CRITICAL (MPI)
845	c$OMP END MASTER
846	#endif
847
848	c$OMP BARRIER
849
850	ENDIF ! using_mpi
851
852
853	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
854	DO l=1,llm
855
856	zdufi2(1:klon,l)=zdufi(1:klon,l)
857	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
858
859	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
860	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
861
862	pdhfi(:,jj_begin,l)=0
863	pdqfi(:,jj_begin,l,:)=0
864	pdufi(:,jj_begin,l)=0
865	pdvfi(:,jj_begin,l)=0
866
867	if (.not. is_south_pole) then
868	pdhfi(:,jj_end,l)=0
869	pdqfi(:,jj_end,l,:)=0
870	pdufi(:,jj_end,l)=0
871	pdvfi(:,jj_end,l)=0
872	endif
873
874	ENDDO
875	c$OMP END DO NOWAIT
876
877	c$OMP MASTER
878	pdpsfi(:,jj_begin)=0
879	if (.not. is_south_pole) then
880	pdpsfi(:,jj_end)=0
881	endif
882	c$OMP END MASTER
883	c-----------------------------------------------------------------------
884	c transformation des tendances physiques en tendances dynamiques:
885	c ---------------------------------------------------------------
886
887	c tendance sur la pression :
888	c -----------------------------------
889	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
890	c
891	c 62. enthalpie potentielle
892	c ---------------------
893
894	kstart=1
895	kend=klon
896
897	if (is_north_pole) kstart=2
898	if (is_south_pole) kend=klon-1
899
900	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
901	DO l=1,llm
902
903	!CDIR ON_ADB(index_i)
904	!CDIR ON_ADB(index_j)
905	!cdir NODEP
906	do ig0=kstart,kend
907	i=index_i(ig0)
908	j=index_j(ig0)
909	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
910	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
911	enddo
912
913	if (is_north_pole) then
914	DO i=1,iip1
915	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
916	enddo
917	endif
918
919	if (is_south_pole) then
920	DO i=1,iip1
921	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
922	ENDDO
923	endif
924	ENDDO
925	c$OMP END DO NOWAIT
926
927	c 62. humidite specifique
928	c ---------------------
929	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
930	! DO iq=1,nqtot
931	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
932	! DO l=1,llm
933	!!!cdir NODEP
934	! do ig0=kstart,kend
935	! i=index_i(ig0)
936	! j=index_j(ig0)
937	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
938	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
939	! enddo
940	!
941	! if (is_north_pole) then
942	! do i=1,iip1
943	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
944	! enddo
945	! endif
946	!
947	! if (is_south_pole) then
948	! do i=1,iip1
949	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
950	! enddo
951	! endif
952	! ENDDO
953	!c$OMP END DO NOWAIT
954	! ENDDO
955
956	c 63. traceurs
957	c ------------
958	C initialisation des tendances
959
960	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
961	DO l=1,llm
962	pdqfi(:,:,l,:)=0.
963	ENDDO
964	c$OMP END DO NOWAIT
965
966	C
967	!cdir NODEP
968	DO iq=1,nqtot
969	iiq=niadv(iq)
970	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
971	DO l=1,llm
972	!CDIR ON_ADB(index_i)
973	!CDIR ON_ADB(index_j)
974	!cdir NODEP
975	DO ig0=kstart,kend
976	i=index_i(ig0)
977	j=index_j(ig0)
978	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
979	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
980	ENDDO
981
982	IF (is_north_pole) then
983	DO i=1,iip1
984	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
985	ENDDO
986	ENDIF
987
988	IF (is_south_pole) then
989	DO i=1,iip1
990	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
991	ENDDO
992	ENDIF
993
994	ENDDO
995	c$OMP END DO NOWAIT
996	ENDDO
997
998	c 65. champ u:
999	c ------------
1000	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1001	DO l=1,llm
1002	!CDIR ON_ADB(index_i)
1003	!CDIR ON_ADB(index_j)
1004	!cdir NODEP
1005	do ig0=kstart,kend
1006	i=index_i(ig0)
1007	j=index_j(ig0)
1008
1009	if (i/=iim) then
1010	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1011	endif
1012
1013	if (i==1) then
1014	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1015	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1016	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1017	endif
1018
1019	enddo
1020
1021	if (is_north_pole) then
1022	DO i=1,iip1
1023	pdufi(i,1,l) = 0.
1024	ENDDO
1025	endif
1026
1027	if (is_south_pole) then
1028	DO i=1,iip1
1029	pdufi(i,jjp1,l) = 0.
1030	ENDDO
1031	endif
1032
1033	ENDDO
1034	c$OMP END DO NOWAIT
1035
1036	c 67. champ v:
1037	c ------------
1038
1039	kstart=1
1040	kend=klon
1041
1042	if (is_north_pole) kstart=2
1043	if (is_south_pole) kend=klon-1-iim
1044
1045	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1046	DO l=1,llm
1047	!CDIR ON_ADB(index_i)
1048	!CDIR ON_ADB(index_j)
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	#else
1113	write(lunout,*)
1114	& "calfis_p: for now can only work with parallel physics"
1115	stop
1116	#endif
1117	! of #ifdef CPP_EARTH
1118	RETURN
1119	END

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

Download in other formats:

Original Format