Context Navigation

calfis_p.F @ 2598

Last change on this file since 2598 was 2597, checked in by Ehouarn Millour, 9 years ago
Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 30.1 KB

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

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

Download in other formats:

Original Format