Context Navigation

calfis_p.F @ 3457

Last change on this file since 3457 was 1393, checked in by Ehouarn Millour, 14 years ago

Some minor fixes for the Newtonian and Shallow Water cases:
-Some initializations were done twice, causing problems with the FFT filter.
-Some CPP_EARTH preprocessing keys should in fact be CPP_IOIPSL.

Also added some useful debugging options to 'debug' mode in arch-PW6_VARGAS.fcm

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 28.0 KB

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

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

Download in other formats:

Original Format