Context Navigation

calfis_loc.F @ 2556

Last change on this file since 2556 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:keywords set to `Id`
File size: 31.5 KB

Line
1	!
2	! $Id: calfis_loc.F 2429 2016-01-27 12:43:09Z jyg $
3	!
4	C
5	C
6	SUBROUTINE calfis_loc(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	c
31	c Auteur : P. Le Van, F. Hourdin
32	c .........
33	USE dimphy
34	USE mod_phys_lmdz_mpi_data, mpi_root_xx=>mpi_master
35	USE mod_phys_lmdz_omp_data, ONLY: klon_omp, klon_omp_begin
36	USE mod_const_mpi, ONLY: COMM_LMDZ
37	USE mod_interface_dyn_phys
38	USE IOPHY
39	#endif
40	#ifdef CPP_PARA
41	USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
42	$ ,jj_begin_dyn=>jj_begin,jj_end_dyn=>jj_end
43	USE Write_Field
44	Use Write_field_p
45	USE Times
46	#endif
47	USE infotrac, ONLY: nqtot, niadv, tname
48	USE control_mod, ONLY: planet_type, nsplit_phys
49	#ifdef CPP_PHYS
50	USE callphysiq_mod, ONLY: call_physiq
51	#endif
52
53	#ifdef CPP_PARA
54	IMPLICIT NONE
55	c=======================================================================
56	c
57	c 1. rearrangement des tableaux et transformation
58	c variables dynamiques > variables physiques
59	c 2. calcul des termes physiques
60	c 3. retransformation des tendances physiques en tendances dynamiques
61	c
62	c remarques:
63	c ----------
64	c
65	c - les vents sont donnes dans la physique par leurs composantes
66	c naturelles.
67	c - la variable thermodynamique de la physique est une variable
68	c intensive : T
69	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
70	c - les deux seules variables dependant de la geometrie necessaires
71	c pour la physique sont la latitude pour le rayonnement et
72	c l'aire de la maille quand on veut integrer une grandeur
73	c horizontalement.
74	c - les points de la physique sont les points scalaires de la
75	c la dynamique; numerotation:
76	c 1 pour le pole nord
77	c (jjm-1)*iim pour l'interieur du domaine
78	c ngridmx pour le pole sud
79	c ---> ngridmx=2+(jjm-1)*iim
80	c
81	c Input :
82	c -------
83	c ecritphy frequence d'ecriture (en jours)de histphy
84	c pucov covariant zonal velocity
85	c pvcov covariant meridional velocity
86	c pteta potential temperature
87	c pps surface pressure
88	c pmasse masse d'air dans chaque maille
89	c pts surface temperature (K)
90	c callrad clef d'appel au rayonnement
91	c
92	c Output :
93	c --------
94	c pdufi tendency for the natural zonal velocity (ms-1)
95	c pdvfi tendency for the natural meridional velocity
96	c pdhfi tendency for the potential temperature
97	c pdtsfi tendency for the surface temperature
98	c
99	c pdtrad radiative tendencies \ both input
100	c pfluxrad radiative fluxes / and output
101	c
102	c=======================================================================
103	c
104	c-----------------------------------------------------------------------
105	c
106	c 0. Declarations :
107	c ------------------
108
109	#include "dimensions.h"
110	#include "paramet.h"
111	#include "temps.h"
112
113	INTEGER ngridmx
114	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
115
116	#include "comconst.h"
117	#include "comvert.h"
118	#include "comgeom2.h"
119	#include "iniprint.h"
120	#ifdef CPP_MPI
121	include 'mpif.h'
122	#endif
123	c Arguments :
124	c -----------
125	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
126	REAL,INTENT(IN):: jD_cur, jH_cur
127	REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
128	REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
129	REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
130	REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
131	REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
132	REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
133	REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential
134
135	REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
136	REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
137	REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
138	REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used
139
140	REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
141	REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
142	REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
143	REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
144
145	! tendencies (in */s) from the physics
146	REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
147	REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
148	REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
149	REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
150	REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)
151
152	#ifdef CPP_PHYS
153	! Ehouarn: for now calfis_p needs some informations from physics to compile
154	c Local variables :
155	c -----------------
156
157	INTEGER i,j,l,ig0,ig,iq,iiq
158	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
159	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
160	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
161	c
162	REAL zrot(iip1,jjb_v:jje_v,llm) ! AdlC May 2014
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	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
375	DO l=1,llm
376	!CDIR ON_ADB(index_i)
377	!CDIR ON_ADB(index_j)
378	do ig0=1,klon
379	i=index_i(ig0)
380	j=index_j(ig0)
381	zphi(ig0,l) = pphi(i,j,l)
382	enddo
383	ENDDO
384	c$OMP END DO NOWAIT
385
386	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
387
388	c$OMP MASTER
389	!CDIR ON_ADB(index_i)
390	!CDIR ON_ADB(index_j)
391	do ig0=1,klon
392	i=index_i(ig0)
393	j=index_j(ig0)
394	zphis(ig0) = pphis(i,j)
395	enddo
396	c$OMP END MASTER
397
398
399	c CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
400
401	c$OMP BARRIER
402
403	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
404	DO l=1,llm
405	DO ig=1,klon
406	zphi(ig,l)=zphi(ig,l)-zphis(ig)
407	ENDDO
408	ENDDO
409	c$OMP END DO NOWAIT
410
411
412	c
413	c 45. champ u:
414	c ------------
415
416	kstart=1
417	kend=klon
418
419	if (is_north_pole_dyn) kstart=2
420	if (is_south_pole_dyn) kend=klon-1
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	!CDIR SPARSE
427	do ig0=kstart,kend
428	i=index_i(ig0)
429	j=index_j(ig0)
430	if (i==1) then
431	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
432	$ + pucov(1,j,l)/cu(1,j) )
433	else
434	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
435	$ + pucov(i,j,l)/cu(i,j) )
436	endif
437	enddo
438	ENDDO
439	c$OMP END DO NOWAIT
440
441	c
442	C Alvaro de la Camara (May 2014)
443	C 46.1 Calcul de la vorticite et passage sur la grille physique
444	C --------------------------------------------------------------
445
446	jjb=jj_begin_dyn-1
447	jje=jj_end_dyn+1
448	if (is_north_pole_dyn) jjb=1
449	if (is_south_pole_dyn) jje=jjm
450
451	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
452
453	DO l=1,llm
454	do i=1,iim
455	do j=jjb,jje
456	zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
457	$ + pucov(i,j+1,l) - pucov(i,j,l))
458	$ / (cu(i,j)+cu(i,j+1))
459	$ / (cv(i+1,j)+cv(i,j)) *4
460	enddo
461	enddo
462	ENDDO
463
464
465	c 46.2champ v:
466	c -----------
467
468	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
469	DO l=1,llm
470	!CDIR ON_ADB(index_i)
471	!CDIR ON_ADB(index_j)
472	DO ig0=kstart,kend
473	i=index_i(ig0)
474	j=index_j(ig0)
475	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
476	$ + pvcov(i,j,l)/cv(i,j) )
477	if (j==1 .OR. j==jjp1) then ! AdlC MAY 2014
478	zrfi(ig0,l) = 0 ! AdlC MAY 2014
479	else
480	if(i==1)then
481	zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
482	$ +zrot(1,j-1,l)+zrot(1,j,l)) ! AdlC MAY 2014
483	else
484	zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
485	$ +zrot(i,j-1,l)+zrot(i,j,l)) ! AdlC MAY 2014
486	endif
487	endif
488
489
490	ENDDO
491	ENDDO
492	c$OMP END DO NOWAIT
493
494	c 47. champs de vents aux pole nord
495	c ------------------------------
496	c U = 1 / pi * integrale [ v * cos(long) * d long ]
497	c V = 1 / pi * integrale [ v * sin(long) * d long ]
498
499	if (is_north_pole_dyn) then
500	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
501	DO l=1,llm
502
503	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
504	DO i=2,iim
505	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
506	ENDDO
507
508	DO i=1,iim
509	zcos(i) = COS(rlonv(i))*z1(i)
510	zsin(i) = SIN(rlonv(i))*z1(i)
511	ENDDO
512
513	zufi(1,l) = SSUM(iim,zcos,1)/pi
514	zvfi(1,l) = SSUM(iim,zsin,1)/pi
515	zrfi(1,l) = 0.
516
517	ENDDO
518	c$OMP END DO NOWAIT
519	endif
520
521
522	c 48. champs de vents aux pole sud:
523	c ---------------------------------
524	c U = 1 / pi * integrale [ v * cos(long) * d long ]
525	c V = 1 / pi * integrale [ v * sin(long) * d long ]
526
527	if (is_south_pole_dyn) then
528	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
529	DO l=1,llm
530
531	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
532	DO i=2,iim
533	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
534	ENDDO
535
536	DO i=1,iim
537	zcos(i) = COS(rlonv(i))*z1(i)
538	zsin(i) = SIN(rlonv(i))*z1(i)
539	ENDDO
540
541	zufi(klon,l) = SSUM(iim,zcos,1)/pi
542	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
543	zrfi(klon,l) = 0.
544	ENDDO
545	c$OMP END DO NOWAIT
546	endif
547
548	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
549	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
550	DO l=1,llm
551	!CDIR ON_ADB(index_i)
552	!CDIR ON_ADB(index_j)
553	do ig0=1,klon
554	i=index_i(ig0)
555	j=index_j(ig0)
556	flxwfi(ig0,l) = flxw(i,j,l)
557	enddo
558	ENDDO
559	c$OMP END DO NOWAIT
560
561	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
562
563	c-----------------------------------------------------------------------
564	c Appel de la physique:
565	c ---------------------
566
567
568	c$OMP BARRIER
569	if (first_omp) then
570	klon=klon_omp
571
572	allocate(zplev_omp(klon,llm+1))
573	allocate(zplay_omp(klon,llm))
574	allocate(zphi_omp(klon,llm))
575	allocate(zphis_omp(klon))
576	allocate(presnivs_omp(llm))
577	allocate(zufi_omp(klon,llm))
578	allocate(zvfi_omp(klon,llm))
579	allocate(zrfi_omp(klon,llm)) ! LG Ari 2014
580	allocate(ztfi_omp(klon,llm))
581	allocate(zqfi_omp(klon,llm,nqtot))
582	allocate(zdufi_omp(klon,llm))
583	allocate(zdvfi_omp(klon,llm))
584	allocate(zdtfi_omp(klon,llm))
585	allocate(zdqfi_omp(klon,llm,nqtot))
586	allocate(zdufic_omp(klon,llm))
587	allocate(zdvfic_omp(klon,llm))
588	allocate(zdtfic_omp(klon,llm))
589	allocate(zdqfic_omp(klon,llm,nqtot))
590	allocate(zdpsrf_omp(klon))
591	allocate(flxwfi_omp(klon,llm))
592	first_omp=.false.
593	endif
594
595
596	klon=klon_omp
597	offset=klon_omp_begin-1
598
599	do l=1,llm+1
600	do i=1,klon
601	zplev_omp(i,l)=zplev(offset+i,l)
602	enddo
603	enddo
604
605	do l=1,llm
606	do i=1,klon
607	zplay_omp(i,l)=zplay(offset+i,l)
608	enddo
609	enddo
610
611	do l=1,llm
612	do i=1,klon
613	zphi_omp(i,l)=zphi(offset+i,l)
614	enddo
615	enddo
616
617	do i=1,klon
618	zphis_omp(i)=zphis(offset+i)
619	enddo
620
621
622	do l=1,llm
623	presnivs_omp(l)=presnivs(l)
624	enddo
625
626	do l=1,llm
627	do i=1,klon
628	zufi_omp(i,l)=zufi(offset+i,l)
629	enddo
630	enddo
631
632	do l=1,llm
633	do i=1,klon
634	zvfi_omp(i,l)=zvfi(offset+i,l)
635	enddo
636	enddo
637
638	do l=1,llm
639	do i=1,klon
640	zrfi_omp(i,l)=zrfi(offset+i,l)
641	enddo
642	enddo
643
644	do l=1,llm
645	do i=1,klon
646	ztfi_omp(i,l)=ztfi(offset+i,l)
647	enddo
648	enddo
649
650	do iq=1,nqtot
651	do l=1,llm
652	do i=1,klon
653	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
654	enddo
655	enddo
656	enddo
657
658	do l=1,llm
659	do i=1,klon
660	zdufi_omp(i,l)=zdufi(offset+i,l)
661	enddo
662	enddo
663
664	do l=1,llm
665	do i=1,klon
666	zdvfi_omp(i,l)=zdvfi(offset+i,l)
667	enddo
668	enddo
669
670	do l=1,llm
671	do i=1,klon
672	zdtfi_omp(i,l)=zdtfi(offset+i,l)
673	enddo
674	enddo
675
676	do iq=1,nqtot
677	do l=1,llm
678	do i=1,klon
679	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
680	enddo
681	enddo
682	enddo
683
684	do i=1,klon
685	zdpsrf_omp(i)=zdpsrf(offset+i)
686	enddo
687
688	do l=1,llm
689	do i=1,klon
690	flxwfi_omp(i,l)=flxwfi(offset+i,l)
691	enddo
692	enddo
693
694	c$OMP BARRIER
695
696
697	!$OMP MASTER
698	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
699	!$OMP END MASTER
700	zdt_split=dtphys/nsplit_phys
701	zdufic_omp(:,:)=0.
702	zdvfic_omp(:,:)=0.
703	zdtfic_omp(:,:)=0.
704	zdqfic_omp(:,:,:)=0.
705
706	#ifdef CPP_PHYS
707	do isplit=1,nsplit_phys
708
709	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
710	debut_split=debut.and.isplit==1
711	lafin_split=lafin.and.isplit==nsplit_phys
712
713	CALL call_physiq(klon,llm,nqtot,tname,
714	& debut_split,lafin_split,
715	& jD_cur,jH_cur_split,zdt_split,
716	& zplev_omp,zplay_omp,
717	& zphi_omp,zphis_omp,
718	& presnivs_omp,
719	& zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
720	& flxwfi_omp,pducov,
721	& zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
722	& zdpsrf_omp)
723
724
725	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
726	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
727	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
728	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
729
730	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
731	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
732	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
733	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
734
735	enddo
736
737	#endif
738	! of #ifdef CPP_PHYS
739
740
741	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
742	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
743	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
744	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
745
746	c$OMP BARRIER
747
748	do l=1,llm+1
749	do i=1,klon
750	zplev(offset+i,l)=zplev_omp(i,l)
751	enddo
752	enddo
753
754	do l=1,llm
755	do i=1,klon
756	zplay(offset+i,l)=zplay_omp(i,l)
757	enddo
758	enddo
759
760	do l=1,llm
761	do i=1,klon
762	zphi(offset+i,l)=zphi_omp(i,l)
763	enddo
764	enddo
765
766
767	do i=1,klon
768	zphis(offset+i)=zphis_omp(i)
769	enddo
770
771
772	do l=1,llm
773	presnivs(l)=presnivs_omp(l)
774	enddo
775
776	do l=1,llm
777	do i=1,klon
778	zufi(offset+i,l)=zufi_omp(i,l)
779	enddo
780	enddo
781
782	do l=1,llm
783	do i=1,klon
784	zvfi(offset+i,l)=zvfi_omp(i,l)
785	enddo
786	enddo
787
788	do l=1,llm
789	do i=1,klon
790	ztfi(offset+i,l)=ztfi_omp(i,l)
791	enddo
792	enddo
793
794	do iq=1,nqtot
795	do l=1,llm
796	do i=1,klon
797	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
798	enddo
799	enddo
800	enddo
801
802	do l=1,llm
803	do i=1,klon
804	zdufi(offset+i,l)=zdufi_omp(i,l)
805	enddo
806	enddo
807
808	do l=1,llm
809	do i=1,klon
810	zdvfi(offset+i,l)=zdvfi_omp(i,l)
811	enddo
812	enddo
813
814	do l=1,llm
815	do i=1,klon
816	zdtfi(offset+i,l)=zdtfi_omp(i,l)
817	enddo
818	enddo
819
820	do iq=1,nqtot
821	do l=1,llm
822	do i=1,klon
823	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
824	enddo
825	enddo
826	enddo
827
828	do i=1,klon
829	zdpsrf(offset+i)=zdpsrf_omp(i)
830	enddo
831
832
833	klon=klon_mpi
834	500 CONTINUE
835	c$OMP BARRIER
836
837	c$OMP MASTER
838	call stop_timer(timer_physic)
839	c$OMP END MASTER
840
841	IF (using_mpi) THEN
842
843	if (MPI_rank>0) then
844
845	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
846	DO l=1,llm
847	du_send(1:iim,l)=zdufi(1:iim,l)
848	dv_send(1:iim,l)=zdvfi(1:iim,l)
849	ENDDO
850	c$OMP END DO NOWAIT
851
852	c$OMP BARRIER
853	#ifdef CPP_MPI
854	c$OMP MASTER
855	!$OMP CRITICAL (MPI)
856	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
857	& COMM_LMDZ,Req(1),ierr)
858	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
859	& COMM_LMDZ,Req(2),ierr)
860	!$OMP END CRITICAL (MPI)
861	c$OMP END MASTER
862	#endif
863	c$OMP BARRIER
864
865	endif
866
867	if (MPI_rank<MPI_Size-1) then
868	c$OMP BARRIER
869	#ifdef CPP_MPI
870	c$OMP MASTER
871	!$OMP CRITICAL (MPI)
872	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
873	& COMM_LMDZ,Req(3),ierr)
874	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
875	& COMM_LMDZ,Req(4),ierr)
876	!$OMP END CRITICAL (MPI)
877	c$OMP END MASTER
878	#endif
879	endif
880
881	c$OMP BARRIER
882
883
884	#ifdef CPP_MPI
885	c$OMP MASTER
886	!$OMP CRITICAL (MPI)
887	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
888	call MPI_WAITALL(4,Req(1),Status,ierr)
889	else if (MPI_rank>0) then
890	call MPI_WAITALL(2,Req(1),Status,ierr)
891	else if (MPI_rank <MPI_Size-1) then
892	call MPI_WAITALL(2,Req(3),Status,ierr)
893	endif
894	!$OMP END CRITICAL (MPI)
895	c$OMP END MASTER
896	#endif
897
898	c$OMP BARRIER
899
900	ENDIF ! using_mpi
901
902
903	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
904	DO l=1,llm
905
906	zdufi2(1:klon,l)=zdufi(1:klon,l)
907	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
908
909	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
910	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
911
912	pdhfi(:,jj_begin,l)=0
913	pdqfi(:,jj_begin,l,:)=0
914	pdufi(:,jj_begin,l)=0
915	pdvfi(:,jj_begin,l)=0
916
917	if (.not. is_south_pole_dyn) then
918	pdhfi(:,jj_end:jj_end+1,l)=0
919	pdqfi(:,jj_end:jj_end+1,l,:)=0
920	pdufi(:,jj_end:jj_end+1,l)=0
921	pdvfi(:,jj_end:jj_end+1,l)=0
922	endif
923
924	ENDDO
925	c$OMP END DO NOWAIT
926
927	c$OMP MASTER
928	pdpsfi(:,jj_begin)=0
929
930	if (.not. is_south_pole_dyn) then
931	pdpsfi(:,jj_end:jj_end+1)=0
932	endif
933	c$OMP END MASTER
934	c-----------------------------------------------------------------------
935	c transformation des tendances physiques en tendances dynamiques:
936	c ---------------------------------------------------------------
937
938	c tendance sur la pression :
939	c -----------------------------------
940	c CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
941
942	c$OMP MASTER
943	kstart=1
944	kend=klon
945
946	if (is_north_pole_dyn) kstart=2
947	if (is_south_pole_dyn) kend=klon-1
948
949	!CDIR ON_ADB(index_i)
950	!CDIR ON_ADB(index_j)
951	!cdir NODEP
952	do ig0=kstart,kend
953	i=index_i(ig0)
954	j=index_j(ig0)
955	pdpsfi(i,j) = zdpsrf(ig0)
956	if (i==1) pdpsfi(iip1,j) = zdpsrf(ig0)
957	enddo
958
959	if (is_north_pole_dyn) then
960	DO i=1,iip1
961	pdpsfi(i,1) = zdpsrf(1)
962	enddo
963	endif
964
965	if (is_south_pole_dyn) then
966	DO i=1,iip1
967	pdpsfi(i,jjp1) = zdpsrf(klon)
968	ENDDO
969	endif
970	c$OMP END MASTER
971	cc$OMP BARRIER
972
973	c
974	c 62. enthalpie potentielle
975	c ---------------------
976
977	kstart=1
978	kend=klon
979
980	if (is_north_pole_dyn) kstart=2
981	if (is_south_pole_dyn) kend=klon-1
982
983	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
984	DO l=1,llm
985
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	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
993	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
994	enddo
995
996	if (is_north_pole_dyn) then
997	DO i=1,iip1
998	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
999	enddo
1000	endif
1001
1002	if (is_south_pole_dyn) then
1003	DO i=1,iip1
1004	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1005	ENDDO
1006	endif
1007	ENDDO
1008	c$OMP END DO NOWAIT
1009
1010	c 62. humidite specifique
1011	c ---------------------
1012	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1013	! DO iq=1,nqtot
1014	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1015	! DO l=1,llm
1016	!!!cdir NODEP
1017	! do ig0=kstart,kend
1018	! i=index_i(ig0)
1019	! j=index_j(ig0)
1020	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1021	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1022	! enddo
1023	!
1024	! if (is_north_pole_dyn) then
1025	! do i=1,iip1
1026	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1027	! enddo
1028	! endif
1029	!
1030	! if (is_south_pole_dyn) then
1031	! do i=1,iip1
1032	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1033	! enddo
1034	! endif
1035	! ENDDO
1036	!c$OMP END DO NOWAIT
1037	! ENDDO
1038
1039	c 63. traceurs
1040	c ------------
1041	C initialisation des tendances
1042
1043	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1044	DO l=1,llm
1045	pdqfi(:,jj_begin:jj_end,l,:)=0.
1046	ENDDO
1047	c$OMP END DO NOWAIT
1048
1049	C
1050	!cdir NODEP
1051	DO iq=1,nqtot
1052	iiq=niadv(iq)
1053	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1054	DO l=1,llm
1055	!CDIR ON_ADB(index_i)
1056	!CDIR ON_ADB(index_j)
1057	!cdir NODEP
1058	DO ig0=kstart,kend
1059	i=index_i(ig0)
1060	j=index_j(ig0)
1061	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1062	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1063	ENDDO
1064
1065	IF (is_north_pole_dyn) then
1066	DO i=1,iip1
1067	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1068	ENDDO
1069	ENDIF
1070
1071	IF (is_south_pole_dyn) then
1072	DO i=1,iip1
1073	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1074	ENDDO
1075	ENDIF
1076
1077	ENDDO
1078	c$OMP END DO NOWAIT
1079	ENDDO
1080
1081	c 65. champ u:
1082	c ------------
1083	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1084	DO l=1,llm
1085	!CDIR ON_ADB(index_i)
1086	!CDIR ON_ADB(index_j)
1087	!cdir NODEP
1088	do ig0=kstart,kend
1089	i=index_i(ig0)
1090	j=index_j(ig0)
1091
1092	if (i/=iim) then
1093	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1094	endif
1095
1096	if (i==1) then
1097	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1098	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1099	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1100	endif
1101
1102	enddo
1103
1104	if (is_north_pole_dyn) then
1105	DO i=1,iip1
1106	pdufi(i,1,l) = 0.
1107	ENDDO
1108	endif
1109
1110	if (is_south_pole_dyn) then
1111	DO i=1,iip1
1112	pdufi(i,jjp1,l) = 0.
1113	ENDDO
1114	endif
1115
1116	ENDDO
1117	c$OMP END DO NOWAIT
1118
1119	c 67. champ v:
1120	c ------------
1121
1122	kstart=1
1123	kend=klon
1124
1125	if (is_north_pole_dyn) kstart=2
1126	if (is_south_pole_dyn) kend=klon-1-iim
1127
1128	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1129	DO l=1,llm
1130	!CDIR ON_ADB(index_i)
1131	!CDIR ON_ADB(index_j)
1132	!cdir NODEP
1133	do ig0=kstart,kend
1134	i=index_i(ig0)
1135	j=index_j(ig0)
1136	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1137	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1138	$ zdvfi2(ig0+iim,l))
1139	$ *cv(i,j)
1140	enddo
1141
1142	ENDDO
1143	c$OMP END DO NOWAIT
1144
1145
1146	c 68. champ v pres des poles:
1147	c ---------------------------
1148	c v = U * cos(long) + V * SIN(long)
1149
1150	if (is_north_pole_dyn) then
1151
1152	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1153	DO l=1,llm
1154
1155	DO i=1,iim
1156	pdvfi(i,1,l)=
1157	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1158
1159	pdvfi(i,1,l)=
1160	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1161	ENDDO
1162
1163	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1164
1165	ENDDO
1166	c$OMP END DO NOWAIT
1167
1168	endif
1169
1170	if (is_south_pole_dyn) then
1171
1172	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1173	DO l=1,llm
1174
1175	DO i=1,iim
1176	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1177	$ +zdvfi(klon,l)*SIN(rlonv(i))
1178
1179	pdvfi(i,jjm,l)=
1180	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1181	ENDDO
1182
1183	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1184
1185	ENDDO
1186	c$OMP END DO NOWAIT
1187
1188	endif
1189	c-----------------------------------------------------------------------
1190
1191	700 CONTINUE
1192
1193	firstcal = .FALSE.
1194
1195	#else
1196	write(lunout,*)
1197	& "calfis_p: for now can only work with parallel physics"
1198	stop
1199	#endif
1200	! of #ifdef CPP_PHYS
1201	#endif
1202	! of #ifdef CPP_PARA
1203	END

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

Download in other formats:

Original Format