Context Navigation

calfis_loc.F @ 4083

Last change on this file since 4083 was 2604, checked in by Ehouarn Millour, 8 years ago
Add Exner function to the call_physiq arguments (not used by the Earth physics) to harmonize physics/dynamics interface. 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:keywords set to `Id`
File size: 32.8 KB

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

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

Download in other formats:

Original Format