Context Navigation

calfis_p.F @ 2031

Last change on this file since 2031 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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 30.3 KB

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

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

Download in other formats:

Original Format