Context Navigation

calfis_p.F @ 2372

Last change on this file since 2372 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

Set things up so that all physics-related initializations are done via iniphysiq.
Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
removed obsolete comgeomphy.h and comgeomphy.F90

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.7 KB

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

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

Download in other formats:

Original Format