Context Navigation

calfis_p.F @ 1249

Last change on this file since 1249 was 1231, checked in by lguez, 15 years ago

These are changes related to ozone in radiative transfer computations:

-- Moved "max(ozonecm, 1e-12)" from procedure "radlwsw" to procedure

"ozonecm".

-- Removed ratio "p_surface / p_reference" in the computation of

"pozon" in "radlwsw".

-- Removed ratio "p_reference / p_surface" in the computation of

"zoz" in "SW_LMDAR4" and "SW_AEROAR4".

-- Removed ratio "mass of air / mass of ozone" in the computation of

"zoz" in "LW_LMDAR4".

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 1231 2009-08-27 12:50:15Z ymeurdesoif $
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, intent(in):: 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