Context Navigation

calfis_loc.F @ 1682

Last change on this file since 1682 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

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

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

Download in other formats:

Original Format