Context Navigation

calfis_loc.F @ 2275

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

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

Download in other formats:

Original Format