Context Navigation

calfis_p.F @ 2225

Last change on this file since 2225 was 2221, checked in by Ehouarn Millour, 10 years ago
Some cleanup: remove (unused) clesph0 from dynamics. 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: 29.5 KB

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

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

Download in other formats:

Original Format