Context Navigation

calfis_loc.F @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1920:1997 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.9 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	cIM diagnostique PVteta, Amip2
222	INTEGER,PARAMETER :: ntetaSTD=3
223	REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
224	REAL PVteta(klon,ntetaSTD)
225
226
227	REAL SSUM
228
229	LOGICAL,SAVE :: firstcal=.true., debut=.true.
230	c$OMP THREADPRIVATE(firstcal,debut)
231
232	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
233	INTEGER :: ierr
234	#ifdef CPP_MPI
235	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
236	#else
237	INTEGER,dimension(1,4) :: Status
238	#endif
239	INTEGER, dimension(4) :: Req
240	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
241	integer :: k,kstart,kend
242	INTEGER :: offset
243
244	LOGICAL tracerdyn
245	c
246	c-----------------------------------------------------------------------
247	c
248	c 1. Initialisations :
249	c --------------------
250	c
251
252	klon=klon_mpi
253
254	PVteta(:,:)=0.
255
256	c
257	IF ( firstcal ) THEN
258	debut = .TRUE.
259	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
260	write(lunout,*) 'STOP dans calfis'
261	write(lunout,*)
262	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
263	write(lunout,*) ' ngridmx jjm iim '
264	write(lunout,*) ngridmx,jjm,iim
265	STOP
266	ENDIF
267	c$OMP MASTER
268	ALLOCATE(zpsrf(klon))
269	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
270	ALLOCATE(zphi(klon,llm),zphis(klon))
271	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
272	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
273	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
274	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
275	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
276	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
277	ALLOCATE(zdpsrf(klon))
278	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
279	ALLOCATE(flxwfi(klon,llm))
280	c$OMP END MASTER
281	c$OMP BARRIER
282	ELSE
283	debut = .FALSE.
284	ENDIF
285
286	c
287	c
288	c-----------------------------------------------------------------------
289	c 40. transformation des variables dynamiques en variables physiques:
290	c ---------------------------------------------------------------
291
292	c 41. pressions au sol (en Pascals)
293	c ----------------------------------
294
295	c$OMP MASTER
296	call start_timer(timer_physic)
297	c$OMP END MASTER
298
299	c$OMP MASTER
300	!CDIR ON_ADB(index_i)
301	!CDIR ON_ADB(index_j)
302	do ig0=1,klon
303	i=index_i(ig0)
304	j=index_j(ig0)
305	zpsrf(ig0)=pps(i,j)
306	enddo
307	c$OMP END MASTER
308
309
310	c 42. pression intercouches :
311	c
312	c -----------------------------------------------------------------
313	c .... zplev definis aux (llm +1) interfaces des couches ....
314	c .... zplay definis aux ( llm ) milieux des couches ....
315	c -----------------------------------------------------------------
316
317	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
318	c
319	unskap = 1./ kappa
320	c
321	c print *,omp_rank,'klon--->',klon
322	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
323	DO l = 1, llmp1
324	!CDIR ON_ADB(index_i)
325	!CDIR ON_ADB(index_j)
326	do ig0=1,klon
327	i=index_i(ig0)
328	j=index_j(ig0)
329	zplev( ig0,l ) = pp(i,j,l)
330	enddo
331	ENDDO
332	c$OMP END DO NOWAIT
333	c
334	c
335
336	c 43. temperature naturelle (en K) et pressions milieux couches .
337	c ---------------------------------------------------------------
338	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
339	DO l=1,llm
340	!CDIR ON_ADB(index_i)
341	!CDIR ON_ADB(index_j)
342	do ig0=1,klon
343	i=index_i(ig0)
344	j=index_j(ig0)
345	pksurcp = ppk(i,j,l) / cpp
346	zplay(ig0,l) = preff * pksurcp ** unskap
347	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
348	enddo
349
350	ENDDO
351	c$OMP END DO NOWAIT
352
353	c 43.bis traceurs
354	c ---------------
355	c
356
357	DO iq=1,nqtot
358	iiq=niadv(iq)
359	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
360	DO l=1,llm
361	!CDIR ON_ADB(index_i)
362	!CDIR ON_ADB(index_j)
363	do ig0=1,klon
364	i=index_i(ig0)
365	j=index_j(ig0)
366	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
367	enddo
368	ENDDO
369	c$OMP END DO NOWAIT
370	ENDDO
371
372
373	c Geopotentiel calcule par rapport a la surface locale:
374	c -----------------------------------------------------
375
376	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
377	DO l=1,llm
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	zphi(ig0,l) = pphi(i,j,l)
384	enddo
385	ENDDO
386	c$OMP END DO NOWAIT
387
388	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
389
390	c$OMP MASTER
391	!CDIR ON_ADB(index_i)
392	!CDIR ON_ADB(index_j)
393	do ig0=1,klon
394	i=index_i(ig0)
395	j=index_j(ig0)
396	zphis(ig0) = pphis(i,j)
397	enddo
398	c$OMP END MASTER
399
400
401	c CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
402
403	c$OMP BARRIER
404
405	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
406	DO l=1,llm
407	DO ig=1,klon
408	zphi(ig,l)=zphi(ig,l)-zphis(ig)
409	ENDDO
410	ENDDO
411	c$OMP END DO NOWAIT
412
413
414	c
415	c 45. champ u:
416	c ------------
417
418	kstart=1
419	kend=klon
420
421	if (is_north_pole) kstart=2
422	if (is_south_pole) kend=klon-1
423
424	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
425	DO l=1,llm
426	!CDIR ON_ADB(index_i)
427	!CDIR ON_ADB(index_j)
428	!CDIR SPARSE
429	do ig0=kstart,kend
430	i=index_i(ig0)
431	j=index_j(ig0)
432	if (i==1) then
433	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
434	$ + pucov(1,j,l)/cu(1,j) )
435	else
436	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
437	$ + pucov(i,j,l)/cu(i,j) )
438	endif
439	enddo
440	ENDDO
441	c$OMP END DO NOWAIT
442
443	c 46.champ v:
444	c -----------
445
446	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
447	DO l=1,llm
448	!CDIR ON_ADB(index_i)
449	!CDIR ON_ADB(index_j)
450	DO ig0=kstart,kend
451	i=index_i(ig0)
452	j=index_j(ig0)
453	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
454	$ + pvcov(i,j,l)/cv(i,j) )
455
456	ENDDO
457	ENDDO
458	c$OMP END DO NOWAIT
459
460	c 47. champs de vents aux pole nord
461	c ------------------------------
462	c U = 1 / pi * integrale [ v * cos(long) * d long ]
463	c V = 1 / pi * integrale [ v * sin(long) * d long ]
464
465	if (is_north_pole) then
466	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
467	DO l=1,llm
468
469	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
470	DO i=2,iim
471	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
472	ENDDO
473
474	DO i=1,iim
475	zcos(i) = COS(rlonv(i))*z1(i)
476	zsin(i) = SIN(rlonv(i))*z1(i)
477	ENDDO
478
479	zufi(1,l) = SSUM(iim,zcos,1)/pi
480	zvfi(1,l) = SSUM(iim,zsin,1)/pi
481
482	ENDDO
483	c$OMP END DO NOWAIT
484	endif
485
486
487	c 48. champs de vents aux pole sud:
488	c ---------------------------------
489	c U = 1 / pi * integrale [ v * cos(long) * d long ]
490	c V = 1 / pi * integrale [ v * sin(long) * d long ]
491
492	if (is_south_pole) then
493	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
494	DO l=1,llm
495
496	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
497	DO i=2,iim
498	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
499	ENDDO
500
501	DO i=1,iim
502	zcos(i) = COS(rlonv(i))*z1(i)
503	zsin(i) = SIN(rlonv(i))*z1(i)
504	ENDDO
505
506	zufi(klon,l) = SSUM(iim,zcos,1)/pi
507	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
508	ENDDO
509	c$OMP END DO NOWAIT
510	endif
511
512
513	IF (is_sequential.and.(planet_type=="earth")) THEN
514	#ifdef CPP_PHYS
515	! PVtheta calls tetalevel, which is in the physics
516	cIM calcul PV a teta=350, 380, 405K
517	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
518	$ ztfi,zplay,zplev,
519	$ ntetaSTD,rtetaSTD,PVteta)
520	c
521	#endif
522	ENDIF
523
524	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
525	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
526	DO l=1,llm
527	!CDIR ON_ADB(index_i)
528	!CDIR ON_ADB(index_j)
529	do ig0=1,klon
530	i=index_i(ig0)
531	j=index_j(ig0)
532	flxwfi(ig0,l) = flxw(i,j,l)
533	enddo
534	ENDDO
535	c$OMP END DO NOWAIT
536
537	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
538
539	c-----------------------------------------------------------------------
540	c Appel de la physique:
541	c ---------------------
542
543
544	c$OMP BARRIER
545	if (first_omp) then
546	klon=klon_omp
547
548	allocate(zplev_omp(klon,llm+1))
549	allocate(zplay_omp(klon,llm))
550	allocate(zphi_omp(klon,llm))
551	allocate(zphis_omp(klon))
552	allocate(presnivs_omp(llm))
553	allocate(zufi_omp(klon,llm))
554	allocate(zvfi_omp(klon,llm))
555	allocate(ztfi_omp(klon,llm))
556	allocate(zqfi_omp(klon,llm,nqtot))
557	allocate(zdufi_omp(klon,llm))
558	allocate(zdvfi_omp(klon,llm))
559	allocate(zdtfi_omp(klon,llm))
560	allocate(zdqfi_omp(klon,llm,nqtot))
561	allocate(zdufic_omp(klon,llm))
562	allocate(zdvfic_omp(klon,llm))
563	allocate(zdtfic_omp(klon,llm))
564	allocate(zdqfic_omp(klon,llm,nqtot))
565	allocate(zdpsrf_omp(klon))
566	allocate(flxwfi_omp(klon,llm))
567	first_omp=.false.
568	endif
569
570
571	klon=klon_omp
572	offset=klon_omp_begin-1
573
574	do l=1,llm+1
575	do i=1,klon
576	zplev_omp(i,l)=zplev(offset+i,l)
577	enddo
578	enddo
579
580	do l=1,llm
581	do i=1,klon
582	zplay_omp(i,l)=zplay(offset+i,l)
583	enddo
584	enddo
585
586	do l=1,llm
587	do i=1,klon
588	zphi_omp(i,l)=zphi(offset+i,l)
589	enddo
590	enddo
591
592	do i=1,klon
593	zphis_omp(i)=zphis(offset+i)
594	enddo
595
596
597	do l=1,llm
598	presnivs_omp(l)=presnivs(l)
599	enddo
600
601	do l=1,llm
602	do i=1,klon
603	zufi_omp(i,l)=zufi(offset+i,l)
604	enddo
605	enddo
606
607	do l=1,llm
608	do i=1,klon
609	zvfi_omp(i,l)=zvfi(offset+i,l)
610	enddo
611	enddo
612
613	do l=1,llm
614	do i=1,klon
615	ztfi_omp(i,l)=ztfi(offset+i,l)
616	enddo
617	enddo
618
619	do iq=1,nqtot
620	do l=1,llm
621	do i=1,klon
622	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
623	enddo
624	enddo
625	enddo
626
627	do l=1,llm
628	do i=1,klon
629	zdufi_omp(i,l)=zdufi(offset+i,l)
630	enddo
631	enddo
632
633	do l=1,llm
634	do i=1,klon
635	zdvfi_omp(i,l)=zdvfi(offset+i,l)
636	enddo
637	enddo
638
639	do l=1,llm
640	do i=1,klon
641	zdtfi_omp(i,l)=zdtfi(offset+i,l)
642	enddo
643	enddo
644
645	do iq=1,nqtot
646	do l=1,llm
647	do i=1,klon
648	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
649	enddo
650	enddo
651	enddo
652
653	do i=1,klon
654	zdpsrf_omp(i)=zdpsrf(offset+i)
655	enddo
656
657	do l=1,llm
658	do i=1,klon
659	flxwfi_omp(i,l)=flxwfi(offset+i,l)
660	enddo
661	enddo
662
663	c$OMP BARRIER
664
665
666	!$OMP MASTER
667	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
668	!$OMP END MASTER
669	zdt_split=dtphys/nsplit_phys
670	zdufic_omp(:,:)=0.
671	zdvfic_omp(:,:)=0.
672	zdtfic_omp(:,:)=0.
673	zdqfic_omp(:,:,:)=0.
674
675	#ifdef CPP_PHYS
676	do isplit=1,nsplit_phys
677
678	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
679	debut_split=debut.and.isplit==1
680	lafin_split=lafin.and.isplit==nsplit_phys
681
682	if (planet_type=="earth") then
683
684	CALL physiq (klon,
685	. llm,
686	. debut_split,
687	. lafin_split,
688	. jD_cur,
689	. jH_cur_split,
690	. zdt_split,
691	. zplev_omp,
692	. zplay_omp,
693	. zphi_omp,
694	. zphis_omp,
695	. presnivs_omp,
696	. clesphy0,
697	. zufi_omp,
698	. zvfi_omp,
699	. ztfi_omp,
700	. zqfi_omp,
701	c#ifdef INCA
702	. flxwfi_omp,
703	c#endif
704	. zdufi_omp,
705	. zdvfi_omp,
706	. zdtfi_omp,
707	. zdqfi_omp,
708	. zdpsrf_omp,
709	cIM diagnostique PVteta, Amip2
710	. pducov,
711	. PVteta)
712
713	else if ( planet_type=="generic" ) then
714
715	CALL physiq (klon, !! ngrid
716	. llm, !! nlayer
717	. nqtot, !! nq
718	. tname, !! tracer names from dynamical core (given in infotrac)
719	. debut_split, !! firstcall
720	. lafin_split, !! lastcall
721	. jD_cur, !! pday. see leapfrog_p
722	. jH_cur_split, !! ptime "fraction of day"
723	. zdt_split, !! ptimestep
724	. zplev_omp, !! pplev
725	. zplay_omp, !! pplay
726	. zphi_omp, !! pphi
727	. zufi_omp, !! pu
728	. zvfi_omp, !! pv
729	. ztfi_omp, !! pt
730	. zqfi_omp, !! pq
731	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
732	. zdufi_omp, !! pdu
733	. zdvfi_omp, !! pdv
734	. zdtfi_omp, !! pdt
735	. zdqfi_omp, !! pdq
736	. zdpsrf_omp, !! pdpsrf
737	. tracerdyn) !! tracerdyn <-- utilite ???
738
739	endif ! of if (planet_type=="earth")
740
741
742	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
743	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
744	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
745	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
746
747	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
748	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
749	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
750	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
751
752	enddo
753
754	#endif
755	! of #ifdef CPP_PHYS
756
757
758	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
759	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
760	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
761	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
762
763	c$OMP BARRIER
764
765	do l=1,llm+1
766	do i=1,klon
767	zplev(offset+i,l)=zplev_omp(i,l)
768	enddo
769	enddo
770
771	do l=1,llm
772	do i=1,klon
773	zplay(offset+i,l)=zplay_omp(i,l)
774	enddo
775	enddo
776
777	do l=1,llm
778	do i=1,klon
779	zphi(offset+i,l)=zphi_omp(i,l)
780	enddo
781	enddo
782
783
784	do i=1,klon
785	zphis(offset+i)=zphis_omp(i)
786	enddo
787
788
789	do l=1,llm
790	presnivs(l)=presnivs_omp(l)
791	enddo
792
793	do l=1,llm
794	do i=1,klon
795	zufi(offset+i,l)=zufi_omp(i,l)
796	enddo
797	enddo
798
799	do l=1,llm
800	do i=1,klon
801	zvfi(offset+i,l)=zvfi_omp(i,l)
802	enddo
803	enddo
804
805	do l=1,llm
806	do i=1,klon
807	ztfi(offset+i,l)=ztfi_omp(i,l)
808	enddo
809	enddo
810
811	do iq=1,nqtot
812	do l=1,llm
813	do i=1,klon
814	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
815	enddo
816	enddo
817	enddo
818
819	do l=1,llm
820	do i=1,klon
821	zdufi(offset+i,l)=zdufi_omp(i,l)
822	enddo
823	enddo
824
825	do l=1,llm
826	do i=1,klon
827	zdvfi(offset+i,l)=zdvfi_omp(i,l)
828	enddo
829	enddo
830
831	do l=1,llm
832	do i=1,klon
833	zdtfi(offset+i,l)=zdtfi_omp(i,l)
834	enddo
835	enddo
836
837	do iq=1,nqtot
838	do l=1,llm
839	do i=1,klon
840	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
841	enddo
842	enddo
843	enddo
844
845	do i=1,klon
846	zdpsrf(offset+i)=zdpsrf_omp(i)
847	enddo
848
849
850	klon=klon_mpi
851	500 CONTINUE
852	c$OMP BARRIER
853
854	c$OMP MASTER
855	call stop_timer(timer_physic)
856	c$OMP END MASTER
857
858	IF (using_mpi) THEN
859
860	if (MPI_rank>0) then
861
862	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
863	DO l=1,llm
864	du_send(1:iim,l)=zdufi(1:iim,l)
865	dv_send(1:iim,l)=zdvfi(1:iim,l)
866	ENDDO
867	c$OMP END DO NOWAIT
868
869	c$OMP BARRIER
870	#ifdef CPP_MPI
871	c$OMP MASTER
872	!$OMP CRITICAL (MPI)
873	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
874	& COMM_LMDZ,Req(1),ierr)
875	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
876	& COMM_LMDZ,Req(2),ierr)
877	!$OMP END CRITICAL (MPI)
878	c$OMP END MASTER
879	#endif
880	c$OMP BARRIER
881
882	endif
883
884	if (MPI_rank<MPI_Size-1) then
885	c$OMP BARRIER
886	#ifdef CPP_MPI
887	c$OMP MASTER
888	!$OMP CRITICAL (MPI)
889	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
890	& COMM_LMDZ,Req(3),ierr)
891	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
892	& COMM_LMDZ,Req(4),ierr)
893	!$OMP END CRITICAL (MPI)
894	c$OMP END MASTER
895	#endif
896	endif
897
898	c$OMP BARRIER
899
900
901	#ifdef CPP_MPI
902	c$OMP MASTER
903	!$OMP CRITICAL (MPI)
904	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
905	call MPI_WAITALL(4,Req(1),Status,ierr)
906	else if (MPI_rank>0) then
907	call MPI_WAITALL(2,Req(1),Status,ierr)
908	else if (MPI_rank <MPI_Size-1) then
909	call MPI_WAITALL(2,Req(3),Status,ierr)
910	endif
911	!$OMP END CRITICAL (MPI)
912	c$OMP END MASTER
913	#endif
914
915	c$OMP BARRIER
916
917	ENDIF ! using_mpi
918
919
920	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
921	DO l=1,llm
922
923	zdufi2(1:klon,l)=zdufi(1:klon,l)
924	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
925
926	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
927	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
928
929	pdhfi(:,jj_begin,l)=0
930	pdqfi(:,jj_begin,l,:)=0
931	pdufi(:,jj_begin,l)=0
932	pdvfi(:,jj_begin,l)=0
933
934	if (.not. is_south_pole) then
935	pdhfi(:,jj_end:jj_end+1,l)=0
936	pdqfi(:,jj_end:jj_end+1,l,:)=0
937	pdufi(:,jj_end:jj_end+1,l)=0
938	pdvfi(:,jj_end:jj_end+1,l)=0
939	endif
940
941	ENDDO
942	c$OMP END DO NOWAIT
943
944	c$OMP MASTER
945	pdpsfi(:,jj_begin)=0
946
947	if (.not. is_south_pole) then
948	pdpsfi(:,jj_end:jj_end+1)=0
949	endif
950	c$OMP END MASTER
951	c-----------------------------------------------------------------------
952	c transformation des tendances physiques en tendances dynamiques:
953	c ---------------------------------------------------------------
954
955	c tendance sur la pression :
956	c -----------------------------------
957	c CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
958
959	c$OMP MASTER
960	kstart=1
961	kend=klon
962
963	if (is_north_pole) kstart=2
964	if (is_south_pole) kend=klon-1
965
966	!CDIR ON_ADB(index_i)
967	!CDIR ON_ADB(index_j)
968	!cdir NODEP
969	do ig0=kstart,kend
970	i=index_i(ig0)
971	j=index_j(ig0)
972	pdpsfi(i,j) = zdpsrf(ig0)
973	if (i==1) pdpsfi(iip1,j) = zdpsrf(ig0)
974	enddo
975
976	if (is_north_pole) then
977	DO i=1,iip1
978	pdpsfi(i,1) = zdpsrf(1)
979	enddo
980	endif
981
982	if (is_south_pole) then
983	DO i=1,iip1
984	pdpsfi(i,jjp1) = zdpsrf(klon)
985	ENDDO
986	endif
987	c$OMP END MASTER
988	cc$OMP BARRIER
989
990	c
991	c 62. enthalpie potentielle
992	c ---------------------
993
994	kstart=1
995	kend=klon
996
997	if (is_north_pole) kstart=2
998	if (is_south_pole) kend=klon-1
999
1000	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1001	DO l=1,llm
1002
1003	!CDIR ON_ADB(index_i)
1004	!CDIR ON_ADB(index_j)
1005	!cdir NODEP
1006	do ig0=kstart,kend
1007	i=index_i(ig0)
1008	j=index_j(ig0)
1009	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1010	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
1011	enddo
1012
1013	if (is_north_pole) then
1014	DO i=1,iip1
1015	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
1016	enddo
1017	endif
1018
1019	if (is_south_pole) then
1020	DO i=1,iip1
1021	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
1022	ENDDO
1023	endif
1024	ENDDO
1025	c$OMP END DO NOWAIT
1026
1027	c 62. humidite specifique
1028	c ---------------------
1029	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
1030	! DO iq=1,nqtot
1031	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1032	! DO l=1,llm
1033	!!!cdir NODEP
1034	! do ig0=kstart,kend
1035	! i=index_i(ig0)
1036	! j=index_j(ig0)
1037	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1038	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1039	! enddo
1040	!
1041	! if (is_north_pole) then
1042	! do i=1,iip1
1043	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1044	! enddo
1045	! endif
1046	!
1047	! if (is_south_pole) then
1048	! do i=1,iip1
1049	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1050	! enddo
1051	! endif
1052	! ENDDO
1053	!c$OMP END DO NOWAIT
1054	! ENDDO
1055
1056	c 63. traceurs
1057	c ------------
1058	C initialisation des tendances
1059
1060	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1061	DO l=1,llm
1062	pdqfi(:,jj_begin:jj_end,l,:)=0.
1063	ENDDO
1064	c$OMP END DO NOWAIT
1065
1066	C
1067	!cdir NODEP
1068	DO iq=1,nqtot
1069	iiq=niadv(iq)
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	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1079	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1080	ENDDO
1081
1082	IF (is_north_pole) then
1083	DO i=1,iip1
1084	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1085	ENDDO
1086	ENDIF
1087
1088	IF (is_south_pole) then
1089	DO i=1,iip1
1090	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1091	ENDDO
1092	ENDIF
1093
1094	ENDDO
1095	c$OMP END DO NOWAIT
1096	ENDDO
1097
1098	c 65. champ u:
1099	c ------------
1100	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1101	DO l=1,llm
1102	!CDIR ON_ADB(index_i)
1103	!CDIR ON_ADB(index_j)
1104	!cdir NODEP
1105	do ig0=kstart,kend
1106	i=index_i(ig0)
1107	j=index_j(ig0)
1108
1109	if (i/=iim) then
1110	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1111	endif
1112
1113	if (i==1) then
1114	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1115	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1116	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1117	endif
1118
1119	enddo
1120
1121	if (is_north_pole) then
1122	DO i=1,iip1
1123	pdufi(i,1,l) = 0.
1124	ENDDO
1125	endif
1126
1127	if (is_south_pole) then
1128	DO i=1,iip1
1129	pdufi(i,jjp1,l) = 0.
1130	ENDDO
1131	endif
1132
1133	ENDDO
1134	c$OMP END DO NOWAIT
1135
1136	c 67. champ v:
1137	c ------------
1138
1139	kstart=1
1140	kend=klon
1141
1142	if (is_north_pole) kstart=2
1143	if (is_south_pole) kend=klon-1-iim
1144
1145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1146	DO l=1,llm
1147	!CDIR ON_ADB(index_i)
1148	!CDIR ON_ADB(index_j)
1149	!cdir NODEP
1150	do ig0=kstart,kend
1151	i=index_i(ig0)
1152	j=index_j(ig0)
1153	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1154	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1155	$ zdvfi2(ig0+iim,l))
1156	$ *cv(i,j)
1157	enddo
1158
1159	ENDDO
1160	c$OMP END DO NOWAIT
1161
1162
1163	c 68. champ v pres des poles:
1164	c ---------------------------
1165	c v = U * cos(long) + V * SIN(long)
1166
1167	if (is_north_pole) then
1168
1169	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1170	DO l=1,llm
1171
1172	DO i=1,iim
1173	pdvfi(i,1,l)=
1174	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1175
1176	pdvfi(i,1,l)=
1177	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1178	ENDDO
1179
1180	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1181
1182	ENDDO
1183	c$OMP END DO NOWAIT
1184
1185	endif
1186
1187	if (is_south_pole) then
1188
1189	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1190	DO l=1,llm
1191
1192	DO i=1,iim
1193	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1194	$ +zdvfi(klon,l)*SIN(rlonv(i))
1195
1196	pdvfi(i,jjm,l)=
1197	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1198	ENDDO
1199
1200	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1201
1202	ENDDO
1203	c$OMP END DO NOWAIT
1204
1205	endif
1206	c-----------------------------------------------------------------------
1207
1208	700 CONTINUE
1209
1210	firstcal = .FALSE.
1211
1212	#else
1213	write(lunout,*)
1214	& "calfis_p: for now can only work with parallel physics"
1215	stop
1216	#endif
1217	! of #ifdef CPP_PHYS
1218	RETURN
1219	END

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

Download in other formats:

Original Format