Context Navigation

calfis_loc.F @ 2242

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

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

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

File size: 31.1 KB

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

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

Download in other formats:

Original Format