Context Navigation

calfis_p.F @ 1227

Last change on this file since 1227 was 1222, checked in by Ehouarn Millour, 15 years ago

Changes and cleanups to enable compiling without physics
and without ioipsl.

IOIPSL related cleanups:

bibio/writehist.F encapsulate the routine (which needs IOIPSL to function)

with #ifdef IOIPSL flag.

dyn3d/abort_gcm.F, dyn3dpar/abort_gcm.F and dyn3dpar/getparam.F90: use ioipsl_getincom module when not compiling with IOIPSL library, in order to always be able to use getin() routine.
removed unused "use IOIPSL" in dyn3dpar/guide_p_mod.F90
calendar related issue: Initialize day_ref and annee_ref in iniacademic.F (i.e. when they are not read from start.nc file)

Earth-specific programs/routines/modules:
create_etat0.F, fluxstokenc.F, limit_netcdf.F, startvar.F
(versions in dyn3d and dyn3dpar)
These routines and modules, which by design and porpose are made to function with
Earth physics are encapsulated with #CPP_EARTH cpp flag.

Earth-specific instructions:

calls to qminimum (specific treatment of first 2 tracers, i.e. water) in dyn3d/caladvtrac.F, dyn3d/integrd.F, dyn3dpar/caladvtrac_p.F, dyn3dpar/integrd_p.F only if (planet_type == 'earth')

Interaction with parallel physics:

routine dyn3dpar/parallel.F90 uses "surface_data" module (which is in the physics ...) to know value of "type_ocean" . Encapsulated that with #ifdef CPP_EARTH and set to a default type_ocean="dummy" otherwise.
So far, only Earth physics are parallelized, so all the interaction between parallel dynamics and parallel physics are encapsulated with #ifdef CCP_EARTH (this way we can run parallel without any physics). The (dyn3dpar) routines which contains such interaction are: bands.F90, gr_dyn_fi_p.F, gr_fi_dyn_p.F, mod_interface_dyn_phys.F90 This should later (when improving dyn/phys interface) be encapsulated with a more general and appropriate #ifdef CPP_PHYS cpp flag.

I checked that these changes do not alter results (on the simple
32x24x11 bench) on Ciclad (seq & mpi), Brodie (seq, mpi & omp) and
Vargas (seq, mpi & omp).

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 28.4 KB

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

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

Download in other formats:

Original Format