Context Navigation

calfis_p.F @ 3485

Last change on this file since 3485 was 2429, checked in by Laurent Fairhead, 9 years ago
Correction on the calculation of the surface of the grid at the poles (problem was introduced in r2222). Due to the different distribution of OMP tasks in the dynamics and the physics, had to introduce 2 new logical variables, is_pole_north_phy and is_pole_south_phy, and so decided to rename the old is_north_pole/is_south_pole to is_north_pole_dyn/is_south_pole_dyn to stay coherent and, hopefully, clear things up a bit. LF
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.0 KB

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

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

Download in other formats:

Original Format