Context Navigation

calfis_loc.F @ 2848

Last change on this file since 2848 was 2037, checked in by lguez, 11 years ago

PVteta computed by PVtheta was not used. Also there were a couple of
problems with PVtheta:

-- PVtheta calls tetalevel in phylmd, and interpolates at

Earth-specific values of potential temperature.

-- PVtheta calls tourabs, and the computation of rot in tourabs should

be modified (it is not correct when there is a zoom).

-- Even when there is no zoom, the computation of rot in tourabs

should probably be modified: directly combine vcov and ucov, as in
tourpot, instead of dividing by cv and cu.

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.3 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	$ clesphy0,
24	$ pdufi,
25	$ pdvfi,
26	$ pdhfi,
27	$ pdqfi,
28	$ pdpsfi)
29	#ifdef CPP_PHYS
30	! If using physics
31	c
32	c Auteur : P. Le Van, F. Hourdin
33	c .........
34	USE dimphy
35	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
36	USE mod_interface_dyn_phys
37	USE IOPHY
38	#endif
39	USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
40	USE Write_Field
41	Use Write_field_p
42	USE Times
43	USE infotrac, ONLY: nqtot, niadv, tname
44	USE control_mod, ONLY: planet_type, nsplit_phys
45
46	IMPLICIT NONE
47	c=======================================================================
48	c
49	c 1. rearrangement des tableaux et transformation
50	c variables dynamiques > variables physiques
51	c 2. calcul des termes physiques
52	c 3. retransformation des tendances physiques en tendances dynamiques
53	c
54	c remarques:
55	c ----------
56	c
57	c - les vents sont donnes dans la physique par leurs composantes
58	c naturelles.
59	c - la variable thermodynamique de la physique est une variable
60	c intensive : T
61	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
62	c - les deux seules variables dependant de la geometrie necessaires
63	c pour la physique sont la latitude pour le rayonnement et
64	c l'aire de la maille quand on veut integrer une grandeur
65	c horizontalement.
66	c - les points de la physique sont les points scalaires de la
67	c la dynamique; numerotation:
68	c 1 pour le pole nord
69	c (jjm-1)*iim pour l'interieur du domaine
70	c ngridmx pour le pole sud
71	c ---> ngridmx=2+(jjm-1)*iim
72	c
73	c Input :
74	c -------
75	c ecritphy frequence d'ecriture (en jours)de histphy
76	c pucov covariant zonal velocity
77	c pvcov covariant meridional velocity
78	c pteta potential temperature
79	c pps surface pressure
80	c pmasse masse d'air dans chaque maille
81	c pts surface temperature (K)
82	c callrad clef d'appel au rayonnement
83	c
84	c Output :
85	c --------
86	c pdufi tendency for the natural zonal velocity (ms-1)
87	c pdvfi tendency for the natural meridional velocity
88	c pdhfi tendency for the potential temperature
89	c pdtsfi tendency for the surface temperature
90	c
91	c pdtrad radiative tendencies \ both input
92	c pfluxrad radiative fluxes / and output
93	c
94	c=======================================================================
95	c
96	c-----------------------------------------------------------------------
97	c
98	c 0. Declarations :
99	c ------------------
100
101	#include "dimensions.h"
102	#include "paramet.h"
103	#include "temps.h"
104
105	INTEGER ngridmx
106	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
107
108	#include "comconst.h"
109	#include "comvert.h"
110	#include "comgeom2.h"
111	#include "iniprint.h"
112	#ifdef CPP_MPI
113	include 'mpif.h'
114	#endif
115	c Arguments :
116	c -----------
117	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
118	REAL,INTENT(IN):: jD_cur, jH_cur
119	REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
120	REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
121	REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
122	REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
123	REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
124	REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
125	REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential
126
127	REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
128	REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
129	REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
130	REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used
131
132	REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
133	REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
134	REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
135	REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm) ! Vertical mass flux on dynamics grid
136
137	! tendencies (in */s) from the physics
138	REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
139	REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
140	REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
141	REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
142	REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)
143
144	INTEGER,PARAMETER :: longcles = 20
145	REAL,INTENT(IN) :: clesphy0( longcles ) ! unused
146
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	. clesphy0,
677	. zufi_omp,
678	. zvfi_omp,
679	. ztfi_omp,
680	. zqfi_omp,
681	c#ifdef INCA
682	. flxwfi_omp,
683	c#endif
684	. zdufi_omp,
685	. zdvfi_omp,
686	. zdtfi_omp,
687	. zdqfi_omp,
688	. zdpsrf_omp,
689	. pducov)
690
691	else if ( planet_type=="generic" ) then
692
693	CALL physiq (klon, !! ngrid
694	. llm, !! nlayer
695	. nqtot, !! nq
696	. tname, !! tracer names from dynamical core (given in infotrac)
697	. debut_split, !! firstcall
698	. lafin_split, !! lastcall
699	. jD_cur, !! pday. see leapfrog_p
700	. jH_cur_split, !! ptime "fraction of day"
701	. zdt_split, !! ptimestep
702	. zplev_omp, !! pplev
703	. zplay_omp, !! pplay
704	. zphi_omp, !! pphi
705	. zufi_omp, !! pu
706	. zvfi_omp, !! pv
707	. ztfi_omp, !! pt
708	. zqfi_omp, !! pq
709	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
710	. zdufi_omp, !! pdu
711	. zdvfi_omp, !! pdv
712	. zdtfi_omp, !! pdt
713	. zdqfi_omp, !! pdq
714	. zdpsrf_omp, !! pdpsrf
715	. tracerdyn) !! tracerdyn <-- utilite ???
716
717	endif ! of if (planet_type=="earth")
718
719
720	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
721	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
722	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
723	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
724
725	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
726	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
727	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
728	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
729
730	enddo
731
732	#endif
733	! of #ifdef CPP_PHYS
734
735
736	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
737	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
738	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
739	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
740
741	c$OMP BARRIER
742
743	do l=1,llm+1
744	do i=1,klon
745	zplev(offset+i,l)=zplev_omp(i,l)
746	enddo
747	enddo
748
749	do l=1,llm
750	do i=1,klon
751	zplay(offset+i,l)=zplay_omp(i,l)
752	enddo
753	enddo
754
755	do l=1,llm
756	do i=1,klon
757	zphi(offset+i,l)=zphi_omp(i,l)
758	enddo
759	enddo
760
761
762	do i=1,klon
763	zphis(offset+i)=zphis_omp(i)
764	enddo
765
766
767	do l=1,llm
768	presnivs(l)=presnivs_omp(l)
769	enddo
770
771	do l=1,llm
772	do i=1,klon
773	zufi(offset+i,l)=zufi_omp(i,l)
774	enddo
775	enddo
776
777	do l=1,llm
778	do i=1,klon
779	zvfi(offset+i,l)=zvfi_omp(i,l)
780	enddo
781	enddo
782
783	do l=1,llm
784	do i=1,klon
785	ztfi(offset+i,l)=ztfi_omp(i,l)
786	enddo
787	enddo
788
789	do iq=1,nqtot
790	do l=1,llm
791	do i=1,klon
792	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
793	enddo
794	enddo
795	enddo
796
797	do l=1,llm
798	do i=1,klon
799	zdufi(offset+i,l)=zdufi_omp(i,l)
800	enddo
801	enddo
802
803	do l=1,llm
804	do i=1,klon
805	zdvfi(offset+i,l)=zdvfi_omp(i,l)
806	enddo
807	enddo
808
809	do l=1,llm
810	do i=1,klon
811	zdtfi(offset+i,l)=zdtfi_omp(i,l)
812	enddo
813	enddo
814
815	do iq=1,nqtot
816	do l=1,llm
817	do i=1,klon
818	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
819	enddo
820	enddo
821	enddo
822
823	do i=1,klon
824	zdpsrf(offset+i)=zdpsrf_omp(i)
825	enddo
826
827
828	klon=klon_mpi
829	500 CONTINUE
830	c$OMP BARRIER
831
832	c$OMP MASTER
833	call stop_timer(timer_physic)
834	c$OMP END MASTER
835
836	IF (using_mpi) THEN
837
838	if (MPI_rank>0) then
839
840	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
841	DO l=1,llm
842	du_send(1:iim,l)=zdufi(1:iim,l)
843	dv_send(1:iim,l)=zdvfi(1:iim,l)
844	ENDDO
845	c$OMP END DO NOWAIT
846
847	c$OMP BARRIER
848	#ifdef CPP_MPI
849	c$OMP MASTER
850	!$OMP CRITICAL (MPI)
851	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
852	& COMM_LMDZ,Req(1),ierr)
853	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
854	& COMM_LMDZ,Req(2),ierr)
855	!$OMP END CRITICAL (MPI)
856	c$OMP END MASTER
857	#endif
858	c$OMP BARRIER
859
860	endif
861
862	if (MPI_rank<MPI_Size-1) then
863	c$OMP BARRIER
864	#ifdef CPP_MPI
865	c$OMP MASTER
866	!$OMP CRITICAL (MPI)
867	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
868	& COMM_LMDZ,Req(3),ierr)
869	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
870	& COMM_LMDZ,Req(4),ierr)
871	!$OMP END CRITICAL (MPI)
872	c$OMP END MASTER
873	#endif
874	endif
875
876	c$OMP BARRIER
877
878
879	#ifdef CPP_MPI
880	c$OMP MASTER
881	!$OMP CRITICAL (MPI)
882	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
883	call MPI_WAITALL(4,Req(1),Status,ierr)
884	else if (MPI_rank>0) then
885	call MPI_WAITALL(2,Req(1),Status,ierr)
886	else if (MPI_rank <MPI_Size-1) then
887	call MPI_WAITALL(2,Req(3),Status,ierr)
888	endif
889	!$OMP END CRITICAL (MPI)
890	c$OMP END MASTER
891	#endif
892
893	c$OMP BARRIER
894
895	ENDIF ! using_mpi
896
897
898	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
899	DO l=1,llm
900
901	zdufi2(1:klon,l)=zdufi(1:klon,l)
902	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
903
904	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
905	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
906
907	pdhfi(:,jj_begin,l)=0
908	pdqfi(:,jj_begin,l,:)=0
909	pdufi(:,jj_begin,l)=0
910	pdvfi(:,jj_begin,l)=0
911
912	if (.not. is_south_pole) then
913	pdhfi(:,jj_end:jj_end+1,l)=0
914	pdqfi(:,jj_end:jj_end+1,l,:)=0
915	pdufi(:,jj_end:jj_end+1,l)=0
916	pdvfi(:,jj_end:jj_end+1,l)=0
917	endif
918
919	ENDDO
920	c$OMP END DO NOWAIT
921
922	c$OMP MASTER
923	pdpsfi(:,jj_begin)=0
924
925	if (.not. is_south_pole) then
926	pdpsfi(:,jj_end:jj_end+1)=0
927	endif
928	c$OMP END MASTER
929	c-----------------------------------------------------------------------
930	c transformation des tendances physiques en tendances dynamiques:
931	c ---------------------------------------------------------------
932
933	c tendance sur la pression :
934	c -----------------------------------
935	c CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
936
937	c$OMP MASTER
938	kstart=1
939	kend=klon
940
941	if (is_north_pole) kstart=2
942	if (is_south_pole) kend=klon-1
943
944	!CDIR ON_ADB(index_i)
945	!CDIR ON_ADB(index_j)
946	!cdir NODEP
947	do ig0=kstart,kend
948	i=index_i(ig0)
949	j=index_j(ig0)
950	pdpsfi(i,j) = zdpsrf(ig0)
951	if (i==1) pdpsfi(iip1,j) = zdpsrf(ig0)
952	enddo
953
954	if (is_north_pole) then
955	DO i=1,iip1
956	pdpsfi(i,1) = zdpsrf(1)
957	enddo
958	endif
959
960	if (is_south_pole) then
961	DO i=1,iip1
962	pdpsfi(i,jjp1) = zdpsrf(klon)
963	ENDDO
964	endif
965	c$OMP END MASTER
966	cc$OMP BARRIER
967
968	c
969	c 62. enthalpie potentielle
970	c ---------------------
971
972	kstart=1
973	kend=klon
974
975	if (is_north_pole) kstart=2
976	if (is_south_pole) kend=klon-1
977
978	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
979	DO l=1,llm
980
981	!CDIR ON_ADB(index_i)
982	!CDIR ON_ADB(index_j)
983	!cdir NODEP
984	do ig0=kstart,kend
985	i=index_i(ig0)
986	j=index_j(ig0)
987	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
988	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
989	enddo
990
991	if (is_north_pole) then
992	DO i=1,iip1
993	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
994	enddo
995	endif
996
997	if (is_south_pole) then
998	DO i=1,iip1
999	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1000	ENDDO
1001	endif
1002	ENDDO
1003	c$OMP END DO NOWAIT
1004
1005	c 62. humidite specifique
1006	c ---------------------
1007	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1008	! DO iq=1,nqtot
1009	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1010	! DO l=1,llm
1011	!!!cdir NODEP
1012	! do ig0=kstart,kend
1013	! i=index_i(ig0)
1014	! j=index_j(ig0)
1015	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1016	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1017	! enddo
1018	!
1019	! if (is_north_pole) then
1020	! do i=1,iip1
1021	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1022	! enddo
1023	! endif
1024	!
1025	! if (is_south_pole) then
1026	! do i=1,iip1
1027	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1028	! enddo
1029	! endif
1030	! ENDDO
1031	!c$OMP END DO NOWAIT
1032	! ENDDO
1033
1034	c 63. traceurs
1035	c ------------
1036	C initialisation des tendances
1037
1038	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1039	DO l=1,llm
1040	pdqfi(:,jj_begin:jj_end,l,:)=0.
1041	ENDDO
1042	c$OMP END DO NOWAIT
1043
1044	C
1045	!cdir NODEP
1046	DO iq=1,nqtot
1047	iiq=niadv(iq)
1048	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1049	DO l=1,llm
1050	!CDIR ON_ADB(index_i)
1051	!CDIR ON_ADB(index_j)
1052	!cdir NODEP
1053	DO ig0=kstart,kend
1054	i=index_i(ig0)
1055	j=index_j(ig0)
1056	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1057	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1058	ENDDO
1059
1060	IF (is_north_pole) then
1061	DO i=1,iip1
1062	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1063	ENDDO
1064	ENDIF
1065
1066	IF (is_south_pole) then
1067	DO i=1,iip1
1068	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1069	ENDDO
1070	ENDIF
1071
1072	ENDDO
1073	c$OMP END DO NOWAIT
1074	ENDDO
1075
1076	c 65. champ u:
1077	c ------------
1078	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1079	DO l=1,llm
1080	!CDIR ON_ADB(index_i)
1081	!CDIR ON_ADB(index_j)
1082	!cdir NODEP
1083	do ig0=kstart,kend
1084	i=index_i(ig0)
1085	j=index_j(ig0)
1086
1087	if (i/=iim) then
1088	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1089	endif
1090
1091	if (i==1) then
1092	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1093	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1094	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1095	endif
1096
1097	enddo
1098
1099	if (is_north_pole) then
1100	DO i=1,iip1
1101	pdufi(i,1,l) = 0.
1102	ENDDO
1103	endif
1104
1105	if (is_south_pole) then
1106	DO i=1,iip1
1107	pdufi(i,jjp1,l) = 0.
1108	ENDDO
1109	endif
1110
1111	ENDDO
1112	c$OMP END DO NOWAIT
1113
1114	c 67. champ v:
1115	c ------------
1116
1117	kstart=1
1118	kend=klon
1119
1120	if (is_north_pole) kstart=2
1121	if (is_south_pole) kend=klon-1-iim
1122
1123	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1124	DO l=1,llm
1125	!CDIR ON_ADB(index_i)
1126	!CDIR ON_ADB(index_j)
1127	!cdir NODEP
1128	do ig0=kstart,kend
1129	i=index_i(ig0)
1130	j=index_j(ig0)
1131	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1132	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1133	$ zdvfi2(ig0+iim,l))
1134	$ *cv(i,j)
1135	enddo
1136
1137	ENDDO
1138	c$OMP END DO NOWAIT
1139
1140
1141	c 68. champ v pres des poles:
1142	c ---------------------------
1143	c v = U * cos(long) + V * SIN(long)
1144
1145	if (is_north_pole) then
1146
1147	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1148	DO l=1,llm
1149
1150	DO i=1,iim
1151	pdvfi(i,1,l)=
1152	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1153
1154	pdvfi(i,1,l)=
1155	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1156	ENDDO
1157
1158	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1159
1160	ENDDO
1161	c$OMP END DO NOWAIT
1162
1163	endif
1164
1165	if (is_south_pole) then
1166
1167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1168	DO l=1,llm
1169
1170	DO i=1,iim
1171	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1172	$ +zdvfi(klon,l)*SIN(rlonv(i))
1173
1174	pdvfi(i,jjm,l)=
1175	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1176	ENDDO
1177
1178	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1179
1180	ENDDO
1181	c$OMP END DO NOWAIT
1182
1183	endif
1184	c-----------------------------------------------------------------------
1185
1186	700 CONTINUE
1187
1188	firstcal = .FALSE.
1189
1190	#else
1191	write(lunout,*)
1192	& "calfis_p: for now can only work with parallel physics"
1193	stop
1194	#endif
1195	! of #ifdef CPP_PHYS
1196	RETURN
1197	END

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

Download in other formats:

Original Format