Context Navigation

calfis_p.F @ 1647

Last change on this file since 1647 was 1617, checked in by Ehouarn Millour, 13 years ago
Minor correction: with previous commits, compilation without physics was no longer possible (because PVtheta calls routine tetalevel, which is in the physics directory). Note also that the levels (hard coded in calfis and calfis_p) given as inputs to PVtheta are Earth-specific. This needs to be cleaned and improved... EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.1 KB

Line
1	!
2	! $Id: calfis_p.F 1617 2012-02-28 11:35:00Z emillour $
3	!
4	C
5	C
6	SUBROUTINE calfis_p(lafin,
7	$ jD_cur, jH_cur,
8	$ pucov,
9	$ pvcov,
10	$ pteta,
11	$ pq,
12	$ pmasse,
13	$ pps,
14	$ pp,
15	$ ppk,
16	$ pphis,
17	$ pphi,
18	$ pducov,
19	$ pdvcov,
20	$ pdteta,
21	$ pdq,
22	$ flxw,
23	$ clesphy0,
24	$ pdufi,
25	$ pdvfi,
26	$ pdhfi,
27	$ pdqfi,
28	$ pdpsfi)
29	#ifdef CPP_PHYS
30	! If using physics
31	USE dimphy
32	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
33	USE mod_interface_dyn_phys
34	USE IOPHY
35	#endif
36	USE parallel, ONLY : omp_chunk, using_mpi
37	USE Write_Field
38	Use Write_field_p
39	USE Times
40	USE infotrac
41	USE control_mod
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 "iniprint.h"
109	#ifdef CPP_MPI
110	include 'mpif.h'
111	#endif
112	c Arguments :
113	c -----------
114	LOGICAL lafin
115	! REAL heure
116	REAL, intent(in):: jD_cur, jH_cur
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	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
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	#ifdef CPP_PHYS
146	! Ehouarn: for now calfis_p needs some informations from physics to compile
147	c Local variables :
148	c -----------------
149
150	INTEGER i,j,l,ig0,ig,iq,iiq
151	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
152	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
153	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
154	c
155	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
156	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
157	c
158	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
159	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
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	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
177	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
178	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
180	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
181	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
182
183	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184	! Introduction du splitting (FH)
185	! Question pour Yann :
186	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
187	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
188	! soit allocatable (plutot par exemple que de passer une dimension
189	! dépendant du process en argument des routines) et que, du coup,
190	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
191	! Tu confirmes ?
192	! J'ai suivi le même principe pour les zdufic_omp
193	! Mais c'est surement bien que tu controles.
194	!
195
196	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
197	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
198	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
199	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
200	REAL jH_cur_split,zdt_split
201	LOGICAL debut_split,lafin_split
202	INTEGER isplit
203	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204
205	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
206	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
207	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
208	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
209	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
210
211	LOGICAL,SAVE :: first_omp=.true.
212	c$OMP THREADPRIVATE(first_omp)
213
214	REAL zsin(iim),zcos(iim),z1(iim)
215	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
216	REAL unskap, pksurcp
217	c
218	cIM diagnostique PVteta, Amip2
219	INTEGER ntetaSTD
220	PARAMETER(ntetaSTD=3)
221	REAL rtetaSTD(ntetaSTD)
222	DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
223	REAL PVteta(klon,ntetaSTD)
224
225
226	REAL SSUM
227
228	LOGICAL firstcal, debut
229	DATA firstcal/.true./
230	SAVE firstcal,debut
231	c$OMP THREADPRIVATE(firstcal,debut)
232
233	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
234	INTEGER :: ierr
235	#ifdef CPP_MPI
236	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
237	#else
238	INTEGER,dimension(1,4) :: Status
239	#endif
240	INTEGER, dimension(4) :: Req
241	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
242	integer :: k,kstart,kend
243	INTEGER :: offset
244	c
245	c-----------------------------------------------------------------------
246	c
247	c 1. Initialisations :
248	c --------------------
249	c
250
251	klon=klon_mpi
252
253	PVteta(:,:)=0.
254
255	c
256	IF ( firstcal ) THEN
257	debut = .TRUE.
258	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
259	write(lunout,*) 'STOP dans calfis'
260	write(lunout,*)
261	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
262	write(lunout,*) ' ngridmx jjm iim '
263	write(lunout,*) 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.and.(planet_type=="earth")) THEN
490	#ifdef CPP_PHYS
491	! PVtheta calls tetalevel, which is in the physics
492	cIM calcul PV a teta=350, 380, 405K
493	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
494	$ ztfi,zplay,zplev,
495	$ ntetaSTD,rtetaSTD,PVteta)
496	c
497	#endif
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	!$OMP MASTER
630	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
631	!$OMP END MASTER
632	zdt_split=dtphys/nsplit_phys
633	zdufic_omp(:,:)=0.
634	zdvfic_omp(:,:)=0.
635	zdtfic_omp(:,:)=0.
636	zdqfic_omp(:,:,:)=0.
637
638	if (planet_type=="earth") then
639	#ifdef CPP_PHYS
640	do isplit=1,nsplit_phys
641
642	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
643	debut_split=debut.and.isplit==1
644	lafin_split=lafin.and.isplit==nsplit_phys
645
646
647	CALL physiq (klon,
648	. llm,
649	. debut_split,
650	. lafin_split,
651	. jD_cur,
652	. jH_cur_split,
653	. zdt_split,
654	. zplev_omp,
655	. zplay_omp,
656	. zphi_omp,
657	. zphis_omp,
658	. presnivs_omp,
659	. clesphy0,
660	. zufi_omp,
661	. zvfi_omp,
662	. ztfi_omp,
663	. zqfi_omp,
664	c#ifdef INCA
665	. flxwfi_omp,
666	c#endif
667	. zdufi_omp,
668	. zdvfi_omp,
669	. zdtfi_omp,
670	. zdqfi_omp,
671	. zdpsrf_omp,
672	cIM diagnostique PVteta, Amip2
673	. pducov,
674	. PVteta)
675
676	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
677	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
678	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
679	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
680
681	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
682	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
683	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
684	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
685
686	enddo
687
688	#endif
689	! of #ifdef CPP_PHYS
690	endif !of if (planet_type=="earth")
691
692	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
693	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
694	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
695	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
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_PHYS
1118	RETURN
1119	END

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

Download in other formats:

Original Format