Context Navigation

calfis_loc.F @ 3831

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

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

Download in other formats:

Original Format