Context Navigation

calfis_p.F @ 832

Last change on this file since 832 was 776, checked in by emillour, 12 years ago
Common dynamics: updates to keep up with LMDZ5 Earth (rev 1649) See file "DOC/chantiers/commit_importants.log" for details. EM
File size: 30.2 KB

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

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

Download in other formats:

Original Format