Context Navigation

calfis_p.F @ 974

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

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

Download in other formats:

Original Format