Context Navigation

calfis_p.F @ 97

Last change on this file since 97 was 97, checked in by slebonnois, 14 years ago
Serie de modifs SL pour homogeneisation des phytitan et phyvenus Ca touche aussi aux liens phy/dyn (surtout a propos de clesphy0), a verifier avec les autres, donc...
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: For now, calfis_p needs Earth 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	USE parallel, ONLY : omp_chunk, using_mpi
39	USE Write_Field
40	Use Write_field_p
41	USE Times
42	USE infotrac
43	USE control_mod
44
45	IMPLICIT NONE
46	c=======================================================================
47	c
48	c 1. rearrangement des tableaux et transformation
49	c variables dynamiques > variables physiques
50	c 2. calcul des termes physiques
51	c 3. retransformation des tendances physiques en tendances dynamiques
52	c
53	c remarques:
54	c ----------
55	c
56	c - les vents sont donnes dans la physique par leurs composantes
57	c naturelles.
58	c - la variable thermodynamique de la physique est une variable
59	c intensive : T
60	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
61	c - les deux seules variables dependant de la geometrie necessaires
62	c pour la physique sont la latitude pour le rayonnement et
63	c l'aire de la maille quand on veut integrer une grandeur
64	c horizontalement.
65	c - les points de la physique sont les points scalaires de la
66	c la dynamique; numerotation:
67	c 1 pour le pole nord
68	c (jjm-1)*iim pour l'interieur du domaine
69	c ngridmx pour le pole sud
70	c ---> ngridmx=2+(jjm-1)*iim
71	c
72	c Input :
73	c -------
74	c ecritphy frequence d'ecriture (en jours)de histphy
75	c pucov covariant zonal velocity
76	c pvcov covariant meridional velocity
77	c pteta potential temperature
78	c pps surface pressure
79	c pmasse masse d'air dans chaque maille
80	c pts surface temperature (K)
81	c callrad clef d'appel au rayonnement
82	c
83	c Output :
84	c --------
85	c pdufi tendency for the natural zonal velocity (ms-1)
86	c pdvfi tendency for the natural meridional velocity
87	c pdhfi tendency for the potential temperature (K/s)
88	c pdtsfi tendency for the surface temperature
89	c
90	c pdtrad radiative tendencies \ both input
91	c pfluxrad radiative fluxes / and output
92	c
93	c=======================================================================
94	c
95	c-----------------------------------------------------------------------
96	c
97	c 0. Declarations :
98	c ------------------
99
100	#include "dimensions.h"
101	#include "paramet.h"
102	#include "temps.h"
103
104	INTEGER ngridmx
105	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
106
107	#include "comconst.h"
108	#include "comvert.h"
109	#include "comgeom2.h"
110	#include "iniprint.h"
111	#ifdef CPP_MPI
112	include 'mpif.h'
113	#endif
114	c Arguments :
115	c -----------
116	LOGICAL lafin
117	! REAL heure
118	REAL, intent(in):: jD_cur, jH_cur
119	REAL pvcov(iip1,jjm,llm)
120	REAL pucov(iip1,jjp1,llm)
121	REAL pteta(iip1,jjp1,llm)
122	REAL pmasse(iip1,jjp1,llm)
123	REAL pq(iip1,jjp1,llm,nqtot)
124	REAL pphis(iip1,jjp1)
125	REAL pphi(iip1,jjp1,llm)
126	c
127	REAL pdvcov(iip1,jjm,llm)
128	REAL pducov(iip1,jjp1,llm)
129	REAL pdteta(iip1,jjp1,llm)
130	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
131	! qui lui meme ne sert a rien dans la routine telle qu'elle est
132	! ecrite, et que j'ai donc commente....
133	REAL pdq(iip1,jjp1,llm,nqtot)
134	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
135	c
136	REAL pps(iip1,jjp1)
137	REAL pp(iip1,jjp1,llmp1)
138	REAL ppk(iip1,jjp1,llm)
139	c
140	REAL pdvfi(iip1,jjm,llm)
141	REAL pdufi(iip1,jjp1,llm)
142	REAL pdhfi(iip1,jjp1,llm)
143	REAL pdqfi(iip1,jjp1,llm,nqtot)
144	REAL pdpsfi(iip1,jjp1)
145
146	#ifdef CPP_PHYS
147	c Local variables :
148	c -----------------
149
150	INTEGER i,j,l,ig0,ig,iq,iiq
151	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
152	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
153	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
154	c
155	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
156	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
157	! ADAPTATION GCM POUR CP(T)
158	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
159	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
160	c
161	! Ces calculs ne servent pas.
162	! Si necessaire, decommenter ces variables et les calculs...
163	! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
164	! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
165	c
166	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
167	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
168	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
169	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
170
171	c
172	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
173	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
174	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
175	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
176	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
177	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
178	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
180	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
181	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
182	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
183	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
184	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
185	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
186	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
187
188	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189	! Introduction du splitting (FH)
190	! Question pour Yann :
191	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
192	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
193	! soit allocatable (plutot par exemple que de passer une dimension
194	! dépendant du process en argument des routines) et que, du coup,
195	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
196	! Tu confirmes ?
197	! J'ai suivi le même principe pour les zdufic_omp
198	! Mais c'est surement bien que tu controles.
199	!
200
201	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
202	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
203	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
204	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
205	REAL jH_cur_split,zdt_split
206	LOGICAL debut_split,lafin_split
207	INTEGER isplit
208	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
209
210	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
211	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
212	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
213	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
214	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
215
216	LOGICAL,SAVE :: first_omp=.true.
217	c$OMP THREADPRIVATE(first_omp)
218
219	REAL zsin(iim),zcos(iim),z1(iim)
220	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
221	REAL unskap, pksurcp
222	c
223	cIM diagnostique PVteta, Amip2
224	INTEGER ntetaSTD
225	PARAMETER(ntetaSTD=3)
226	REAL rtetaSTD(ntetaSTD)
227	DATA rtetaSTD/350., 380., 405./
228	REAL PVteta(klon,ntetaSTD)
229
230
231	REAL SSUM
232
233	LOGICAL firstcal, debut
234	DATA firstcal/.true./
235	SAVE firstcal,debut
236	c$OMP THREADPRIVATE(firstcal,debut)
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	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(ngridmx*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) THEN
515	c
516	cIM calcul PV a teta=350, 380, 405K
517	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
518	$ ztfi,zplay,zplev,
519	$ ntetaSTD,rtetaSTD,PVteta)
520	c
521	ENDIF
522
523	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
524	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
525
526	c-----------------------------------------------------------------------
527	c Appel de la physique:
528	c ---------------------
529
530	! Appel de la physique: pose probleme quand on tourne
531	! SANS physique, car physiq.F est dans le repertoire phy[]...
532	! Il faut une cle CPP_PHYS
533
534	! Le fait que les arguments de physiq soient differents selon les planetes
535	! ne pose pas de probleme a priori.
536
537
538	c$OMP BARRIER
539	if (first_omp) then
540	klon=klon_omp
541
542	allocate(zplev_omp(klon,llm+1))
543	allocate(zplay_omp(klon,llm))
544	allocate(zphi_omp(klon,llm))
545	allocate(zphis_omp(klon))
546	allocate(presnivs_omp(llm))
547	allocate(zufi_omp(klon,llm))
548	allocate(zvfi_omp(klon,llm))
549	allocate(ztfi_omp(klon,llm))
550	allocate(zqfi_omp(klon,llm,nqtot))
551	allocate(zdufi_omp(klon,llm))
552	allocate(zdvfi_omp(klon,llm))
553	allocate(zdtfi_omp(klon,llm))
554	allocate(zdqfi_omp(klon,llm,nqtot))
555	allocate(zdufic_omp(klon,llm))
556	allocate(zdvfic_omp(klon,llm))
557	allocate(zdtfic_omp(klon,llm))
558	allocate(zdqfic_omp(klon,llm,nqtot))
559	allocate(zdpsrf_omp(klon))
560	allocate(flxwfi_omp(klon,llm))
561	first_omp=.false.
562	endif
563
564
565	klon=klon_omp
566	offset=klon_omp_begin-1
567
568	do l=1,llm+1
569	do i=1,klon
570	zplev_omp(i,l)=zplev(offset+i,l)
571	enddo
572	enddo
573
574	do l=1,llm
575	do i=1,klon
576	zplay_omp(i,l)=zplay(offset+i,l)
577	enddo
578	enddo
579
580	do l=1,llm
581	do i=1,klon
582	zphi_omp(i,l)=zphi(offset+i,l)
583	enddo
584	enddo
585
586	do i=1,klon
587	zphis_omp(i)=zphis(offset+i)
588	enddo
589
590
591	do l=1,llm
592	presnivs_omp(l)=presnivs(l)
593	enddo
594
595	do l=1,llm
596	do i=1,klon
597	zufi_omp(i,l)=zufi(offset+i,l)
598	enddo
599	enddo
600
601	do l=1,llm
602	do i=1,klon
603	zvfi_omp(i,l)=zvfi(offset+i,l)
604	enddo
605	enddo
606
607	do l=1,llm
608	do i=1,klon
609	ztfi_omp(i,l)=ztfi(offset+i,l)
610	enddo
611	enddo
612
613	do iq=1,nqtot
614	do l=1,llm
615	do i=1,klon
616	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
617	enddo
618	enddo
619	enddo
620
621	do l=1,llm
622	do i=1,klon
623	zdufi_omp(i,l)=zdufi(offset+i,l)
624	enddo
625	enddo
626
627	do l=1,llm
628	do i=1,klon
629	zdvfi_omp(i,l)=zdvfi(offset+i,l)
630	enddo
631	enddo
632
633	do l=1,llm
634	do i=1,klon
635	zdtfi_omp(i,l)=zdtfi(offset+i,l)
636	enddo
637	enddo
638
639	do iq=1,nqtot
640	do l=1,llm
641	do i=1,klon
642	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
643	enddo
644	enddo
645	enddo
646
647	do i=1,klon
648	zdpsrf_omp(i)=zdpsrf(offset+i)
649	enddo
650
651	do l=1,llm
652	do i=1,klon
653	flxwfi_omp(i,l)=flxwfi(offset+i,l)
654	enddo
655	enddo
656
657	c$OMP BARRIER
658
659	!$OMP MASTER
660	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
661	!$OMP END MASTER
662	zdt_split=dtphys/nsplit_phys
663	zdufic_omp(:,:)=0.
664	zdvfic_omp(:,:)=0.
665	zdtfic_omp(:,:)=0.
666	zdqfic_omp(:,:,:)=0.
667
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	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
745	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
746	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
747	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
748
749	c$OMP BARRIER
750
751	do l=1,llm+1
752	do i=1,klon
753	zplev(offset+i,l)=zplev_omp(i,l)
754	enddo
755	enddo
756
757	do l=1,llm
758	do i=1,klon
759	zplay(offset+i,l)=zplay_omp(i,l)
760	enddo
761	enddo
762
763	do l=1,llm
764	do i=1,klon
765	zphi(offset+i,l)=zphi_omp(i,l)
766	enddo
767	enddo
768
769
770	do i=1,klon
771	zphis(offset+i)=zphis_omp(i)
772	enddo
773
774
775	do l=1,llm
776	presnivs(l)=presnivs_omp(l)
777	enddo
778
779	do l=1,llm
780	do i=1,klon
781	zufi(offset+i,l)=zufi_omp(i,l)
782	enddo
783	enddo
784
785	do l=1,llm
786	do i=1,klon
787	zvfi(offset+i,l)=zvfi_omp(i,l)
788	enddo
789	enddo
790
791	do l=1,llm
792	do i=1,klon
793	ztfi(offset+i,l)=ztfi_omp(i,l)
794	enddo
795	enddo
796
797	do iq=1,nqtot
798	do l=1,llm
799	do i=1,klon
800	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
801	enddo
802	enddo
803	enddo
804
805	do l=1,llm
806	do i=1,klon
807	zdufi(offset+i,l)=zdufi_omp(i,l)
808	enddo
809	enddo
810
811	do l=1,llm
812	do i=1,klon
813	zdvfi(offset+i,l)=zdvfi_omp(i,l)
814	enddo
815	enddo
816
817	do l=1,llm
818	do i=1,klon
819	zdtfi(offset+i,l)=zdtfi_omp(i,l)
820	enddo
821	enddo
822
823	do iq=1,nqtot
824	do l=1,llm
825	do i=1,klon
826	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
827	enddo
828	enddo
829	enddo
830
831	do i=1,klon
832	zdpsrf(offset+i)=zdpsrf_omp(i)
833	enddo
834
835
836	klon=klon_mpi
837	500 CONTINUE
838	c$OMP BARRIER
839
840	c$OMP MASTER
841	call stop_timer(timer_physic)
842	c$OMP END MASTER
843
844	IF (using_mpi) THEN
845
846	if (MPI_rank>0) then
847
848	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
849	DO l=1,llm
850	du_send(1:iim,l)=zdufi(1:iim,l)
851	dv_send(1:iim,l)=zdvfi(1:iim,l)
852	ENDDO
853	c$OMP END DO NOWAIT
854
855	c$OMP BARRIER
856	#ifdef CPP_MPI
857	c$OMP MASTER
858	!$OMP CRITICAL (MPI)
859	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
860	& COMM_LMDZ,Req(1),ierr)
861	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
862	& COMM_LMDZ,Req(2),ierr)
863	!$OMP END CRITICAL (MPI)
864	c$OMP END MASTER
865	#endif
866	c$OMP BARRIER
867
868	endif
869
870	if (MPI_rank<MPI_Size-1) then
871	c$OMP BARRIER
872	#ifdef CPP_MPI
873	c$OMP MASTER
874	!$OMP CRITICAL (MPI)
875	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
876	& COMM_LMDZ,Req(3),ierr)
877	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
878	& COMM_LMDZ,Req(4),ierr)
879	!$OMP END CRITICAL (MPI)
880	c$OMP END MASTER
881	#endif
882	endif
883
884	c$OMP BARRIER
885
886
887	#ifdef CPP_MPI
888	c$OMP MASTER
889	!$OMP CRITICAL (MPI)
890	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
891	call MPI_WAITALL(4,Req(1),Status,ierr)
892	else if (MPI_rank>0) then
893	call MPI_WAITALL(2,Req(1),Status,ierr)
894	else if (MPI_rank <MPI_Size-1) then
895	call MPI_WAITALL(2,Req(3),Status,ierr)
896	endif
897	!$OMP END CRITICAL (MPI)
898	c$OMP END MASTER
899	#endif
900
901	c$OMP BARRIER
902
903	ENDIF ! using_mpi
904
905
906	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
907	DO l=1,llm
908
909	zdufi2(1:klon,l)=zdufi(1:klon,l)
910	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
911
912	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
913	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
914
915	pdhfi(:,jj_begin,l)=0
916	pdqfi(:,jj_begin,l,:)=0
917	pdufi(:,jj_begin,l)=0
918	pdvfi(:,jj_begin,l)=0
919
920	if (.not. is_south_pole) then
921	pdhfi(:,jj_end,l)=0
922	pdqfi(:,jj_end,l,:)=0
923	pdufi(:,jj_end,l)=0
924	pdvfi(:,jj_end,l)=0
925	endif
926
927	ENDDO
928	c$OMP END DO NOWAIT
929
930	c$OMP MASTER
931	pdpsfi(:,jj_begin)=0
932	if (.not. is_south_pole) then
933	pdpsfi(:,jj_end)=0
934	endif
935	c$OMP END MASTER
936	c-----------------------------------------------------------------------
937	c transformation des tendances physiques en tendances dynamiques:
938	c ---------------------------------------------------------------
939
940	c tendance sur la pression :
941	c -----------------------------------
942	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
943	c
944	c 62. enthalpie potentielle
945	c ---------------------
946
947	kstart=1
948	kend=klon
949
950	if (is_north_pole) kstart=2
951	if (is_south_pole) kend=klon-1
952
953	! ADAPTATION GCM POUR CP(T)
954	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
955	DO l=1,llm
956	ztfi(1:klon,l)=ztfi(1:klon,l)+zdtfi(1:klon,l)*dtphys
957	ENDDO
958	!$OMP END DO
959	call t2tpot_p(ngridmx,llm,ztfi,zteta,zpk)
960
961
962	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
963	DO l=1,llm
964	!CDIR ON_ADB(index_i)
965	!CDIR ON_ADB(index_j)
966	!cdir NODEP
967	do ig0=kstart,kend
968	i=index_i(ig0)
969	j=index_j(ig0)
970	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
971	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
972	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
973	if (i==1) then
974	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
975	endif
976	enddo
977
978	if (is_north_pole) then
979	DO i=1,iip1
980	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
981	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
982	enddo
983	endif
984
985	if (is_south_pole) then
986	DO i=1,iip1
987	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
988	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
989	ENDDO
990	endif
991	ENDDO
992	c$OMP END DO NOWAIT
993
994	c 62. humidite specifique
995	c ---------------------
996	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
997	! DO iq=1,nqtot
998	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
999	! DO l=1,llm
1000	!!!cdir NODEP
1001	! do ig0=kstart,kend
1002	! i=index_i(ig0)
1003	! j=index_j(ig0)
1004	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1005	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1006	! enddo
1007	!
1008	! if (is_north_pole) then
1009	! do i=1,iip1
1010	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1011	! enddo
1012	! endif
1013	!
1014	! if (is_south_pole) then
1015	! do i=1,iip1
1016	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1017	! enddo
1018	! endif
1019	! ENDDO
1020	!c$OMP END DO NOWAIT
1021	! ENDDO
1022
1023	c 63. traceurs (tous en intensifs)
1024	c ------------
1025	C initialisation des tendances
1026
1027	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1028	DO l=1,llm
1029	pdqfi(:,:,l,:)=0.
1030	ENDDO
1031	c$OMP END DO NOWAIT
1032
1033	C
1034	!cdir NODEP
1035	DO iq=1,nqtot
1036	iiq=niadv(iq)
1037	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1038	DO l=1,llm
1039	!CDIR ON_ADB(index_i)
1040	!CDIR ON_ADB(index_j)
1041	!cdir NODEP
1042	DO ig0=kstart,kend
1043	i=index_i(ig0)
1044	j=index_j(ig0)
1045	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1046	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1047	ENDDO
1048
1049	IF (is_north_pole) then
1050	DO i=1,iip1
1051	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1052	ENDDO
1053	ENDIF
1054
1055	IF (is_south_pole) then
1056	DO i=1,iip1
1057	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1058	ENDDO
1059	ENDIF
1060
1061	ENDDO
1062	c$OMP END DO NOWAIT
1063	ENDDO
1064
1065	c 65. champ u:
1066	c ------------
1067	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1068	DO l=1,llm
1069	!CDIR ON_ADB(index_i)
1070	!CDIR ON_ADB(index_j)
1071	!cdir NODEP
1072	do ig0=kstart,kend
1073	i=index_i(ig0)
1074	j=index_j(ig0)
1075
1076	if (i/=iim) then
1077	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1078	endif
1079
1080	if (i==1) then
1081	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1082	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1083	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1084	endif
1085
1086	enddo
1087
1088	if (is_north_pole) then
1089	DO i=1,iip1
1090	pdufi(i,1,l) = 0.
1091	ENDDO
1092	endif
1093
1094	if (is_south_pole) then
1095	DO i=1,iip1
1096	pdufi(i,jjp1,l) = 0.
1097	ENDDO
1098	endif
1099
1100	ENDDO
1101	c$OMP END DO NOWAIT
1102
1103	c 67. champ v:
1104	c ------------
1105
1106	kstart=1
1107	kend=klon
1108
1109	if (is_north_pole) kstart=2
1110	if (is_south_pole) kend=klon-1-iim
1111
1112	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1113	DO l=1,llm
1114	!CDIR ON_ADB(index_i)
1115	!CDIR ON_ADB(index_j)
1116	!cdir NODEP
1117	do ig0=kstart,kend
1118	i=index_i(ig0)
1119	j=index_j(ig0)
1120	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1121	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1122	$ zdvfi2(ig0+iim,l))
1123	$ *cv(i,j)
1124	enddo
1125
1126	ENDDO
1127	c$OMP END DO NOWAIT
1128
1129
1130	c 68. champ v pres des poles:
1131	c ---------------------------
1132	c v = U * cos(long) + V * SIN(long)
1133
1134	if (is_north_pole) then
1135
1136	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1137	DO l=1,llm
1138
1139	DO i=1,iim
1140	pdvfi(i,1,l)=
1141	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1142
1143	pdvfi(i,1,l)=
1144	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1145	ENDDO
1146
1147	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1148
1149	ENDDO
1150	c$OMP END DO NOWAIT
1151
1152	endif
1153
1154	if (is_south_pole) then
1155
1156	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1157	DO l=1,llm
1158
1159	DO i=1,iim
1160	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1161	$ +zdvfi(klon,l)*SIN(rlonv(i))
1162
1163	pdvfi(i,jjm,l)=
1164	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1165	ENDDO
1166
1167	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1168
1169	ENDDO
1170	c$OMP END DO NOWAIT
1171
1172	endif
1173	c-----------------------------------------------------------------------
1174
1175	700 CONTINUE
1176
1177	firstcal = .FALSE.
1178
1179	#else
1180	write(lunout,*)
1181	& "calfis_p: for now can only work with parallel physics"
1182	stop
1183	#endif
1184	! of #ifdef CPP_PHYS
1185	RETURN
1186	END

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

Download in other formats:

Original Format