Context Navigation

calfis_loc.F @ 2422

Last change on this file since 2422 was 2418, checked in by Ehouarn Millour, 9 years ago

Improving the physics/dynamics interface:

added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
turned "physiq.F90" into module "physiq_mod.F90" for better control of "physiq" arguments. Extracted embeded "gr_fi_ecrit" as self-standing routine (but note that this routine actually only works in serial mode).

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

File size: 31.4 KB

Line
1	!
2	! $Id$
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) kstart=2
420	if (is_south_pole) 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) jjb=1
449	if (is_south_pole) 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) 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) 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) 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) 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) kstart=2
947	if (is_south_pole) 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) then
960	DO i=1,iip1
961	pdpsfi(i,1) = zdpsrf(1)
962	enddo
963	endif
964
965	if (is_south_pole) 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) kstart=2
981	if (is_south_pole) 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) 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) 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) 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) 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) 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) 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) then
1105	DO i=1,iip1
1106	pdufi(i,1,l) = 0.
1107	ENDDO
1108	endif
1109
1110	if (is_south_pole) 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) kstart=2
1126	if (is_south_pole) 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) 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) 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