Context Navigation

calfis_loc.F @ 1997

Last change on this file since 1997 was 1987, checked in by Ehouarn Millour, 11 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
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