Context Navigation

calfis_loc.F @ 2372

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

More on physics/dynamics separation and cleanup:

Set things up so that all physics-related initializations are done via iniphysiq.
Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
removed obsolete comgeomphy.h and comgeomphy.F90

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

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

Download in other formats:

Original Format