Context Navigation

calfis_loc.F @ 1669

Last change on this file since 1669 was 1669, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1668

File size: 29.7 KB

Line
1	!
2	! $Id$
3	!
4	C
5	C
6	SUBROUTINE calfis_loc(lafin,
7	$ jD_cur, jH_cur,
8	$ pucov,
9	$ pvcov,
10	$ pteta,
11	$ pq,
12	$ pmasse,
13	$ pps,
14	$ pp,
15	$ ppk,
16	$ pphis,
17	$ pphi,
18	$ pducov,
19	$ pdvcov,
20	$ pdteta,
21	$ pdq,
22	$ flxw,
23	$ clesphy0,
24	$ pdufi,
25	$ pdvfi,
26	$ pdhfi,
27	$ pdqfi,
28	$ pdpsfi)
29	#ifdef CPP_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 parallel, ONLY : omp_chunk, using_mpi,jjb_u,jje_u,jjb_v,jje_v
37	USE mod_interface_dyn_phys
38	USE Write_Field
39	Use Write_field_p
40	USE Times
41	USE IOPHY
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
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
119	REAL pvcov(iip1,jjb_v:jje_v,llm)
120	REAL pucov(iip1,jjb_u:jje_u,llm)
121	REAL pteta(iip1,jjb_u:jje_u,llm)
122	REAL pmasse(iip1,jjb_u:jje_u,llm)
123	REAL pq(iip1,jjb_u:jje_u,llm,nqtot)
124	REAL pphis(iip1,jjb_u:jje_u)
125	REAL pphi(iip1,jjb_u:jje_u,llm)
126	c
127	REAL pdvcov(iip1,jjb_v:jje_v,llm)
128	REAL pducov(iip1,jjb_u:jje_u,llm)
129	REAL pdteta(iip1,jjb_u:jje_u,llm)
130	REAL pdq(iip1,jjb_u:jje_u,llm,nqtot)
131	c
132	REAL pps(iip1,jjb_u:jje_u)
133	REAL pp(iip1,jjb_u:jje_u,llmp1)
134	REAL ppk(iip1,jjb_u:jje_u,llm)
135	c
136	REAL pdvfi(iip1,jjb_v:jje_v,llm)
137	REAL pdufi(iip1,jjb_u:jje_u,llm)
138	REAL pdhfi(iip1,jjb_u:jje_u,llm)
139	REAL pdqfi(iip1,jjb_u:jje_u,llm,nqtot)
140	REAL pdpsfi(iip1,jjb_u:jje_u)
141
142	INTEGER longcles
143	PARAMETER ( longcles = 20 )
144	REAL clesphy0( longcles )
145
146
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	c
158	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
159	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
160	c
161	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
162	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
163	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
164	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
165
166	c
167	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
168	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
169	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
170	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
171	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
172	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
173	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
174	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
175	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
176	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
177	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
178	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
180	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
181	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
182
183	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184	! Introduction du splitting (FH)
185	! Question pour Yann :
186	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
187	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
188	! soit allocatable (plutot par exemple que de passer une dimension
189	! dépendant du process en argument des routines) et que, du coup,
190	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
191	! Tu confirmes ?
192	! J'ai suivi le même principe pour les zdufic_omp
193	! Mais c'est surement bien que tu controles.
194	!
195
196	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
197	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
198	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
199	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
200	REAL jH_cur_split,zdt_split
201	LOGICAL debut_split,lafin_split
202	INTEGER isplit
203	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204
205	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
206	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
207	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
208	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
209	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
210
211	LOGICAL,SAVE :: first_omp=.true.
212	c$OMP THREADPRIVATE(first_omp)
213
214	REAL zsin(iim),zcos(iim),z1(iim)
215	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
216	REAL unskap, pksurcp
217	c
218	cIM diagnostique PVteta, Amip2
219	INTEGER ntetaSTD
220	PARAMETER(ntetaSTD=3)
221	REAL rtetaSTD(ntetaSTD)
222	DATA rtetaSTD/350., 380., 405./
223	REAL PVteta(klon,ntetaSTD)
224
225	REAL flxw(iip1,jjb_u:jje_u,llm) ! Flux de masse verticale sur la grille dynamique
226
227	REAL SSUM
228
229	LOGICAL firstcal, debut
230	DATA firstcal/.true./
231	SAVE firstcal,debut
232	c$OMP THREADPRIVATE(firstcal,debut)
233	REAL, intent(in):: jD_cur, jH_cur
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	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
378	DO l=1,llm
379	!CDIR ON_ADB(index_i)
380	!CDIR ON_ADB(index_j)
381	do ig0=1,klon
382	i=index_i(ig0)
383	j=index_j(ig0)
384	zphi(ig0,l) = pphi(i,j,l)
385	enddo
386	ENDDO
387	c$OMP END DO NOWAIT
388
389	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
390
391	c$OMP MASTER
392	!CDIR ON_ADB(index_i)
393	!CDIR ON_ADB(index_j)
394	do ig0=1,klon
395	i=index_i(ig0)
396	j=index_j(ig0)
397	zphis(ig0) = pphis(i,j)
398	enddo
399	c$OMP END MASTER
400
401
402	c 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	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
525	DO l=1,llm
526	!CDIR ON_ADB(index_i)
527	!CDIR ON_ADB(index_j)
528	do ig0=1,klon
529	i=index_i(ig0)
530	j=index_j(ig0)
531	flxwfi(ig0,l) = flxw(i,j,l)
532	enddo
533	ENDDO
534	c$OMP END DO NOWAIT
535
536	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
537
538	c-----------------------------------------------------------------------
539	c Appel de la physique:
540	c ---------------------
541
542
543	c$OMP BARRIER
544	if (first_omp) then
545	klon=klon_omp
546
547	allocate(zplev_omp(klon,llm+1))
548	allocate(zplay_omp(klon,llm))
549	allocate(zphi_omp(klon,llm))
550	allocate(zphis_omp(klon))
551	allocate(presnivs_omp(llm))
552	allocate(zufi_omp(klon,llm))
553	allocate(zvfi_omp(klon,llm))
554	allocate(ztfi_omp(klon,llm))
555	allocate(zqfi_omp(klon,llm,nqtot))
556	allocate(zdufi_omp(klon,llm))
557	allocate(zdvfi_omp(klon,llm))
558	allocate(zdtfi_omp(klon,llm))
559	allocate(zdqfi_omp(klon,llm,nqtot))
560	allocate(zdufic_omp(klon,llm))
561	allocate(zdvfic_omp(klon,llm))
562	allocate(zdtfic_omp(klon,llm))
563	allocate(zdqfic_omp(klon,llm,nqtot))
564	allocate(zdpsrf_omp(klon))
565	allocate(flxwfi_omp(klon,llm))
566	first_omp=.false.
567	endif
568
569
570	klon=klon_omp
571	offset=klon_omp_begin-1
572
573	do l=1,llm+1
574	do i=1,klon
575	zplev_omp(i,l)=zplev(offset+i,l)
576	enddo
577	enddo
578
579	do l=1,llm
580	do i=1,klon
581	zplay_omp(i,l)=zplay(offset+i,l)
582	enddo
583	enddo
584
585	do l=1,llm
586	do i=1,klon
587	zphi_omp(i,l)=zphi(offset+i,l)
588	enddo
589	enddo
590
591	do i=1,klon
592	zphis_omp(i)=zphis(offset+i)
593	enddo
594
595
596	do l=1,llm
597	presnivs_omp(l)=presnivs(l)
598	enddo
599
600	do l=1,llm
601	do i=1,klon
602	zufi_omp(i,l)=zufi(offset+i,l)
603	enddo
604	enddo
605
606	do l=1,llm
607	do i=1,klon
608	zvfi_omp(i,l)=zvfi(offset+i,l)
609	enddo
610	enddo
611
612	do l=1,llm
613	do i=1,klon
614	ztfi_omp(i,l)=ztfi(offset+i,l)
615	enddo
616	enddo
617
618	do iq=1,nqtot
619	do l=1,llm
620	do i=1,klon
621	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
622	enddo
623	enddo
624	enddo
625
626	do l=1,llm
627	do i=1,klon
628	zdufi_omp(i,l)=zdufi(offset+i,l)
629	enddo
630	enddo
631
632	do l=1,llm
633	do i=1,klon
634	zdvfi_omp(i,l)=zdvfi(offset+i,l)
635	enddo
636	enddo
637
638	do l=1,llm
639	do i=1,klon
640	zdtfi_omp(i,l)=zdtfi(offset+i,l)
641	enddo
642	enddo
643
644	do iq=1,nqtot
645	do l=1,llm
646	do i=1,klon
647	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
648	enddo
649	enddo
650	enddo
651
652	do i=1,klon
653	zdpsrf_omp(i)=zdpsrf(offset+i)
654	enddo
655
656	do l=1,llm
657	do i=1,klon
658	flxwfi_omp(i,l)=flxwfi(offset+i,l)
659	enddo
660	enddo
661
662	c$OMP BARRIER
663
664	if (planet_type=="earth") then
665	#ifdef CPP_EARTH
666
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	do isplit=1,nsplit_phys
678
679	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
680	debut_split=debut.and.isplit==1
681	lafin_split=lafin.and.isplit==nsplit_phys
682
683
684	CALL physiq (klon,
685	. llm,
686	. debut_split,
687	. lafin_split,
688	. jD_cur,
689	. jH_cur_split,
690	. zdt_split,
691	. zplev_omp,
692	. zplay_omp,
693	. zphi_omp,
694	. zphis_omp,
695	. presnivs_omp,
696	. clesphy0,
697	. zufi_omp,
698	. zvfi_omp,
699	. ztfi_omp,
700	. zqfi_omp,
701	c#ifdef INCA
702	. flxwfi_omp,
703	c#endif
704	. zdufi_omp,
705	. zdvfi_omp,
706	. zdtfi_omp,
707	. zdqfi_omp,
708	. zdpsrf_omp,
709	cIM diagnostique PVteta, Amip2
710	. pducov,
711	. PVteta)
712
713	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
714	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
715	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
716	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
717
718	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
719	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
720	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
721	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
722
723	enddo
724
725	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
726	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
727	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
728	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
729
730	#endif
731	endif !of if (planet_type=="earth")
732	c$OMP BARRIER
733
734	do l=1,llm+1
735	do i=1,klon
736	zplev(offset+i,l)=zplev_omp(i,l)
737	enddo
738	enddo
739
740	do l=1,llm
741	do i=1,klon
742	zplay(offset+i,l)=zplay_omp(i,l)
743	enddo
744	enddo
745
746	do l=1,llm
747	do i=1,klon
748	zphi(offset+i,l)=zphi_omp(i,l)
749	enddo
750	enddo
751
752
753	do i=1,klon
754	zphis(offset+i)=zphis_omp(i)
755	enddo
756
757
758	do l=1,llm
759	presnivs(l)=presnivs_omp(l)
760	enddo
761
762	do l=1,llm
763	do i=1,klon
764	zufi(offset+i,l)=zufi_omp(i,l)
765	enddo
766	enddo
767
768	do l=1,llm
769	do i=1,klon
770	zvfi(offset+i,l)=zvfi_omp(i,l)
771	enddo
772	enddo
773
774	do l=1,llm
775	do i=1,klon
776	ztfi(offset+i,l)=ztfi_omp(i,l)
777	enddo
778	enddo
779
780	do iq=1,nqtot
781	do l=1,llm
782	do i=1,klon
783	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
784	enddo
785	enddo
786	enddo
787
788	do l=1,llm
789	do i=1,klon
790	zdufi(offset+i,l)=zdufi_omp(i,l)
791	enddo
792	enddo
793
794	do l=1,llm
795	do i=1,klon
796	zdvfi(offset+i,l)=zdvfi_omp(i,l)
797	enddo
798	enddo
799
800	do l=1,llm
801	do i=1,klon
802	zdtfi(offset+i,l)=zdtfi_omp(i,l)
803	enddo
804	enddo
805
806	do iq=1,nqtot
807	do l=1,llm
808	do i=1,klon
809	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
810	enddo
811	enddo
812	enddo
813
814	do i=1,klon
815	zdpsrf(offset+i)=zdpsrf_omp(i)
816	enddo
817
818
819	klon=klon_mpi
820	500 CONTINUE
821	c$OMP BARRIER
822
823	c$OMP MASTER
824	call stop_timer(timer_physic)
825	c$OMP END MASTER
826
827	IF (using_mpi) THEN
828
829	if (MPI_rank>0) then
830
831	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
832	DO l=1,llm
833	du_send(1:iim,l)=zdufi(1:iim,l)
834	dv_send(1:iim,l)=zdvfi(1:iim,l)
835	ENDDO
836	c$OMP END DO NOWAIT
837
838	c$OMP BARRIER
839	#ifdef CPP_MPI
840	c$OMP MASTER
841	!$OMP CRITICAL (MPI)
842	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
843	& COMM_LMDZ,Req(1),ierr)
844	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
845	& COMM_LMDZ,Req(2),ierr)
846	!$OMP END CRITICAL (MPI)
847	c$OMP END MASTER
848	#endif
849	c$OMP BARRIER
850
851	endif
852
853	if (MPI_rank<MPI_Size-1) then
854	c$OMP BARRIER
855	#ifdef CPP_MPI
856	c$OMP MASTER
857	!$OMP CRITICAL (MPI)
858	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
859	& COMM_LMDZ,Req(3),ierr)
860	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
861	& COMM_LMDZ,Req(4),ierr)
862	!$OMP END CRITICAL (MPI)
863	c$OMP END MASTER
864	#endif
865	endif
866
867	c$OMP BARRIER
868
869
870	#ifdef CPP_MPI
871	c$OMP MASTER
872	!$OMP CRITICAL (MPI)
873	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
874	call MPI_WAITALL(4,Req(1),Status,ierr)
875	else if (MPI_rank>0) then
876	call MPI_WAITALL(2,Req(1),Status,ierr)
877	else if (MPI_rank <MPI_Size-1) then
878	call MPI_WAITALL(2,Req(3),Status,ierr)
879	endif
880	!$OMP END CRITICAL (MPI)
881	c$OMP END MASTER
882	#endif
883
884	c$OMP BARRIER
885
886	ENDIF ! using_mpi
887
888
889	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
890	DO l=1,llm
891
892	zdufi2(1:klon,l)=zdufi(1:klon,l)
893	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
894
895	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
896	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
897
898	pdhfi(:,jj_begin,l)=0
899	pdqfi(:,jj_begin,l,:)=0
900	pdufi(:,jj_begin,l)=0
901	pdvfi(:,jj_begin,l)=0
902
903	if (.not. is_south_pole) then
904	pdhfi(:,jj_end:jj_end+1,l)=0
905	pdqfi(:,jj_end:jj_end+1,l,:)=0
906	pdufi(:,jj_end:jj_end+1,l)=0
907	pdvfi(:,jj_end:jj_end+1,l)=0
908	endif
909
910	ENDDO
911	c$OMP END DO NOWAIT
912
913	c$OMP MASTER
914	pdpsfi(:,jj_begin)=0
915
916	if (.not. is_south_pole) then
917	pdpsfi(:,jj_end:jj_end+1)=0
918	endif
919	c$OMP END MASTER
920	c-----------------------------------------------------------------------
921	c transformation des tendances physiques en tendances dynamiques:
922	c ---------------------------------------------------------------
923
924	c tendance sur la pression :
925	c -----------------------------------
926	c CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
927
928	c$OMP MASTER
929	kstart=1
930	kend=klon
931
932	if (is_north_pole) kstart=2
933	if (is_south_pole) kend=klon-1
934
935	!CDIR ON_ADB(index_i)
936	!CDIR ON_ADB(index_j)
937	!cdir NODEP
938	do ig0=kstart,kend
939	i=index_i(ig0)
940	j=index_j(ig0)
941	pdpsfi(i,j) = zdpsrf(ig0)
942	if (i==1) pdpsfi(iip1,j) = zdpsrf(ig0)
943	enddo
944
945	if (is_north_pole) then
946	DO i=1,iip1
947	pdpsfi(i,1) = zdpsrf(1)
948	enddo
949	endif
950
951	if (is_south_pole) then
952	DO i=1,iip1
953	pdpsfi(i,jjp1) = zdpsrf(klon)
954	ENDDO
955	endif
956	c$OMP END MASTER
957	cc$OMP BARRIER
958
959	c
960	c 62. enthalpie potentielle
961	c ---------------------
962
963	kstart=1
964	kend=klon
965
966	if (is_north_pole) kstart=2
967	if (is_south_pole) kend=klon-1
968
969	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
970	DO l=1,llm
971
972	!CDIR ON_ADB(index_i)
973	!CDIR ON_ADB(index_j)
974	!cdir NODEP
975	do ig0=kstart,kend
976	i=index_i(ig0)
977	j=index_j(ig0)
978	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
979	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
980	enddo
981
982	if (is_north_pole) then
983	DO i=1,iip1
984	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
985	enddo
986	endif
987
988	if (is_south_pole) then
989	DO i=1,iip1
990	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
991	ENDDO
992	endif
993	ENDDO
994	c$OMP END DO NOWAIT
995
996	c 62. humidite specifique
997	c ---------------------
998	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
999	! DO iq=1,nqtot
1000	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1001	! DO l=1,llm
1002	!!!cdir NODEP
1003	! do ig0=kstart,kend
1004	! i=index_i(ig0)
1005	! j=index_j(ig0)
1006	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
1007	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
1008	! enddo
1009	!
1010	! if (is_north_pole) then
1011	! do i=1,iip1
1012	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
1013	! enddo
1014	! endif
1015	!
1016	! if (is_south_pole) then
1017	! do i=1,iip1
1018	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
1019	! enddo
1020	! endif
1021	! ENDDO
1022	!c$OMP END DO NOWAIT
1023	! ENDDO
1024
1025	c 63. traceurs
1026	c ------------
1027	C initialisation des tendances
1028
1029	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1030	DO l=1,llm
1031	pdqfi(:,jj_begin:jj_end,l,:)=0.
1032	ENDDO
1033	c$OMP END DO NOWAIT
1034
1035	C
1036	!cdir NODEP
1037	DO iq=1,nqtot
1038	iiq=niadv(iq)
1039	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1040	DO l=1,llm
1041	!CDIR ON_ADB(index_i)
1042	!CDIR ON_ADB(index_j)
1043	!cdir NODEP
1044	DO ig0=kstart,kend
1045	i=index_i(ig0)
1046	j=index_j(ig0)
1047	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1048	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1049	ENDDO
1050
1051	IF (is_north_pole) then
1052	DO i=1,iip1
1053	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1054	ENDDO
1055	ENDIF
1056
1057	IF (is_south_pole) then
1058	DO i=1,iip1
1059	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1060	ENDDO
1061	ENDIF
1062
1063	ENDDO
1064	c$OMP END DO NOWAIT
1065	ENDDO
1066
1067	c 65. champ u:
1068	c ------------
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
1078	if (i/=iim) then
1079	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1080	endif
1081
1082	if (i==1) then
1083	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1084	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1085	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1086	endif
1087
1088	enddo
1089
1090	if (is_north_pole) then
1091	DO i=1,iip1
1092	pdufi(i,1,l) = 0.
1093	ENDDO
1094	endif
1095
1096	if (is_south_pole) then
1097	DO i=1,iip1
1098	pdufi(i,jjp1,l) = 0.
1099	ENDDO
1100	endif
1101
1102	ENDDO
1103	c$OMP END DO NOWAIT
1104
1105	c 67. champ v:
1106	c ------------
1107
1108	kstart=1
1109	kend=klon
1110
1111	if (is_north_pole) kstart=2
1112	if (is_south_pole) kend=klon-1-iim
1113
1114	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1115	DO l=1,llm
1116	!CDIR ON_ADB(index_i)
1117	!CDIR ON_ADB(index_j)
1118	!cdir NODEP
1119	do ig0=kstart,kend
1120	i=index_i(ig0)
1121	j=index_j(ig0)
1122	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1123	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1124	$ zdvfi2(ig0+iim,l))
1125	$ *cv(i,j)
1126	enddo
1127
1128	ENDDO
1129	c$OMP END DO NOWAIT
1130
1131
1132	c 68. champ v pres des poles:
1133	c ---------------------------
1134	c v = U * cos(long) + V * SIN(long)
1135
1136	if (is_north_pole) then
1137
1138	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1139	DO l=1,llm
1140
1141	DO i=1,iim
1142	pdvfi(i,1,l)=
1143	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1144
1145	pdvfi(i,1,l)=
1146	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1147	ENDDO
1148
1149	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1150
1151	ENDDO
1152	c$OMP END DO NOWAIT
1153
1154	endif
1155
1156	if (is_south_pole) then
1157
1158	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1159	DO l=1,llm
1160
1161	DO i=1,iim
1162	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1163	$ +zdvfi(klon,l)*SIN(rlonv(i))
1164
1165	pdvfi(i,jjm,l)=
1166	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1167	ENDDO
1168
1169	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1170
1171	ENDDO
1172	c$OMP END DO NOWAIT
1173
1174	endif
1175	c-----------------------------------------------------------------------
1176
1177	700 CONTINUE
1178
1179	firstcal = .FALSE.
1180
1181	#else
1182	write(,) "calfis_p: for now can only work with parallel physics"
1183	write(lunout,*)
1184	& "calfis_p: for now can only work with parallel physics"
1185	stop
1186	#endif
1187	! of #ifdef CPP_EARTH
1188	RETURN
1189	END

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

Download in other formats:

Original Format