Context Navigation

calfis_p.F @ 5440

Last change on this file since 5440 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
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 29.7 KB

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

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

Download in other formats:

Original Format