Context Navigation

calfis_p.F @ 5452

Last change on this file since 5452 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

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

File size: 29.3 KB

Line
1	!
2	! $Id: calfis_p.F 1676 2012-11-08 20:03:30Z aborella $
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_PHYS
30	! If using physics
31	USE dimphy
32	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
33	USE mod_interface_dyn_phys
34	USE IOPHY
35	#endif
36	USE parallel, ONLY : omp_chunk, using_mpi
37	USE Write_Field
38	Use Write_field_p
39	USE Times
40	USE infotrac
41	USE control_mod
42
43	IMPLICIT NONE
44	c=======================================================================
45	c
46	c 1. rearrangement des tableaux et transformation
47	c variables dynamiques > variables physiques
48	c 2. calcul des termes physiques
49	c 3. retransformation des tendances physiques en tendances dynamiques
50	c
51	c remarques:
52	c ----------
53	c
54	c - les vents sont donnes dans la physique par leurs composantes
55	c naturelles.
56	c - la variable thermodynamique de la physique est une variable
57	c intensive : T
58	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
59	c - les deux seules variables dependant de la geometrie necessaires
60	c pour la physique sont la latitude pour le rayonnement et
61	c l'aire de la maille quand on veut integrer une grandeur
62	c horizontalement.
63	c - les points de la physique sont les points scalaires de la
64	c la dynamique; numerotation:
65	c 1 pour le pole nord
66	c (jjm-1)*iim pour l'interieur du domaine
67	c ngridmx pour le pole sud
68	c ---> ngridmx=2+(jjm-1)*iim
69	c
70	c Input :
71	c -------
72	c ecritphy frequence d'ecriture (en jours)de histphy
73	c pucov covariant zonal velocity
74	c pvcov covariant meridional velocity
75	c pteta potential temperature
76	c pps surface pressure
77	c pmasse masse d'air dans chaque maille
78	c pts surface temperature (K)
79	c callrad clef d'appel au rayonnement
80	c
81	c Output :
82	c --------
83	c pdufi tendency for the natural zonal velocity (ms-1)
84	c pdvfi tendency for the natural meridional velocity
85	c pdhfi tendency for the potential temperature
86	c pdtsfi tendency for the surface temperature
87	c
88	c pdtrad radiative tendencies \ both input
89	c pfluxrad radiative fluxes / and output
90	c
91	c=======================================================================
92	c
93	c-----------------------------------------------------------------------
94	c
95	c 0. Declarations :
96	c ------------------
97
98	#include "dimensions.h"
99	#include "paramet.h"
100	#include "temps.h"
101
102	INTEGER ngridmx
103	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
104
105	#include "comconst.h"
106	#include "comvert.h"
107	#include "comgeom2.h"
108	#include "iniprint.h"
109	#ifdef CPP_MPI
110	include 'mpif.h'
111	#endif
112	c Arguments :
113	c -----------
114	LOGICAL lafin
115	! REAL heure
116	REAL, intent(in):: jD_cur, jH_cur
117	REAL pvcov(iip1,jjm,llm)
118	REAL pucov(iip1,jjp1,llm)
119	REAL pteta(iip1,jjp1,llm)
120	REAL pmasse(iip1,jjp1,llm)
121	REAL pq(iip1,jjp1,llm,nqtot)
122	REAL pphis(iip1,jjp1)
123	REAL pphi(iip1,jjp1,llm)
124	c
125	REAL pdvcov(iip1,jjm,llm)
126	REAL pducov(iip1,jjp1,llm)
127	REAL pdteta(iip1,jjp1,llm)
128	REAL pdq(iip1,jjp1,llm,nqtot)
129	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
130	c
131	REAL pps(iip1,jjp1)
132	REAL pp(iip1,jjp1,llmp1)
133	REAL ppk(iip1,jjp1,llm)
134	c
135	REAL pdvfi(iip1,jjm,llm)
136	REAL pdufi(iip1,jjp1,llm)
137	REAL pdhfi(iip1,jjp1,llm)
138	REAL pdqfi(iip1,jjp1,llm,nqtot)
139	REAL pdpsfi(iip1,jjp1)
140
141	INTEGER longcles
142	PARAMETER ( longcles = 20 )
143	REAL clesphy0( longcles )
144
145	#ifdef CPP_PHYS
146	! Ehouarn: for now calfis_p needs some informations from physics to compile
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./ ! Earth-specific values, beware !!
223	REAL PVteta(klon,ntetaSTD)
224
225
226	REAL SSUM
227
228	LOGICAL firstcal, debut
229	DATA firstcal/.true./
230	SAVE firstcal,debut
231	c$OMP THREADPRIVATE(firstcal,debut)
232
233	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
234	INTEGER :: ierr
235	#ifdef CPP_MPI
236	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
237	#else
238	INTEGER,dimension(1,4) :: Status
239	#endif
240	INTEGER, dimension(4) :: Req
241	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
242	integer :: k,kstart,kend
243	INTEGER :: offset
244
245	LOGICAL tracerdyn
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.and.(planet_type=="earth")) THEN
492	#ifdef CPP_PHYS
493	! PVtheta calls tetalevel, which is in the physics
494	cIM calcul PV a teta=350, 380, 405K
495	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
496	$ ztfi,zplay,zplev,
497	$ ntetaSTD,rtetaSTD,PVteta)
498	c
499	#endif
500	ENDIF
501
502	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
503	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
504
505	c-----------------------------------------------------------------------
506	c Appel de la physique:
507	c ---------------------
508
509
510	c$OMP BARRIER
511	if (first_omp) then
512	klon=klon_omp
513
514	allocate(zplev_omp(klon,llm+1))
515	allocate(zplay_omp(klon,llm))
516	allocate(zphi_omp(klon,llm))
517	allocate(zphis_omp(klon))
518	allocate(presnivs_omp(llm))
519	allocate(zufi_omp(klon,llm))
520	allocate(zvfi_omp(klon,llm))
521	allocate(ztfi_omp(klon,llm))
522	allocate(zqfi_omp(klon,llm,nqtot))
523	allocate(zdufi_omp(klon,llm))
524	allocate(zdvfi_omp(klon,llm))
525	allocate(zdtfi_omp(klon,llm))
526	allocate(zdqfi_omp(klon,llm,nqtot))
527	allocate(zdufic_omp(klon,llm))
528	allocate(zdvfic_omp(klon,llm))
529	allocate(zdtfic_omp(klon,llm))
530	allocate(zdqfic_omp(klon,llm,nqtot))
531	allocate(zdpsrf_omp(klon))
532	allocate(flxwfi_omp(klon,llm))
533	first_omp=.false.
534	endif
535
536
537	klon=klon_omp
538	offset=klon_omp_begin-1
539
540	do l=1,llm+1
541	do i=1,klon
542	zplev_omp(i,l)=zplev(offset+i,l)
543	enddo
544	enddo
545
546	do l=1,llm
547	do i=1,klon
548	zplay_omp(i,l)=zplay(offset+i,l)
549	enddo
550	enddo
551
552	do l=1,llm
553	do i=1,klon
554	zphi_omp(i,l)=zphi(offset+i,l)
555	enddo
556	enddo
557
558	do i=1,klon
559	zphis_omp(i)=zphis(offset+i)
560	enddo
561
562
563	do l=1,llm
564	presnivs_omp(l)=presnivs(l)
565	enddo
566
567	do l=1,llm
568	do i=1,klon
569	zufi_omp(i,l)=zufi(offset+i,l)
570	enddo
571	enddo
572
573	do l=1,llm
574	do i=1,klon
575	zvfi_omp(i,l)=zvfi(offset+i,l)
576	enddo
577	enddo
578
579	do l=1,llm
580	do i=1,klon
581	ztfi_omp(i,l)=ztfi(offset+i,l)
582	enddo
583	enddo
584
585	do iq=1,nqtot
586	do l=1,llm
587	do i=1,klon
588	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
589	enddo
590	enddo
591	enddo
592
593	do l=1,llm
594	do i=1,klon
595	zdufi_omp(i,l)=zdufi(offset+i,l)
596	enddo
597	enddo
598
599	do l=1,llm
600	do i=1,klon
601	zdvfi_omp(i,l)=zdvfi(offset+i,l)
602	enddo
603	enddo
604
605	do l=1,llm
606	do i=1,klon
607	zdtfi_omp(i,l)=zdtfi(offset+i,l)
608	enddo
609	enddo
610
611	do iq=1,nqtot
612	do l=1,llm
613	do i=1,klon
614	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
615	enddo
616	enddo
617	enddo
618
619	do i=1,klon
620	zdpsrf_omp(i)=zdpsrf(offset+i)
621	enddo
622
623	do l=1,llm
624	do i=1,klon
625	flxwfi_omp(i,l)=flxwfi(offset+i,l)
626	enddo
627	enddo
628
629	c$OMP BARRIER
630
631	!$OMP MASTER
632	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
633	!$OMP END MASTER
634	zdt_split=dtphys/nsplit_phys
635	zdufic_omp(:,:)=0.
636	zdvfic_omp(:,:)=0.
637	zdtfic_omp(:,:)=0.
638	zdqfic_omp(:,:,:)=0.
639
640	#ifdef CPP_PHYS
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	if (planet_type=="earth") then
648
649	CALL physiq (klon,
650	. llm,
651	. debut_split,
652	. lafin_split,
653	. jD_cur,
654	. jH_cur_split,
655	. zdt_split,
656	. zplev_omp,
657	. zplay_omp,
658	. zphi_omp,
659	. zphis_omp,
660	. presnivs_omp,
661	. clesphy0,
662	. zufi_omp,
663	. zvfi_omp,
664	. ztfi_omp,
665	. zqfi_omp,
666	c#ifdef INCA
667	. flxwfi_omp,
668	c#endif
669	. zdufi_omp,
670	. zdvfi_omp,
671	. zdtfi_omp,
672	. zdqfi_omp,
673	. zdpsrf_omp,
674	cIM diagnostique PVteta, Amip2
675	. pducov,
676	. PVteta)
677
678	else if ( planet_type=="generic" ) then
679
680	CALL physiq (klon, !! ngrid
681	. llm, !! nlayer
682	. nqtot, !! nq
683	. tname, !! tracer names from dynamical core (given in infotrac)
684	. debut_split, !! firstcall
685	. lafin_split, !! lastcall
686	. jD_cur, !! pday. see leapfrog_p
687	. jH_cur_split, !! ptime "fraction of day"
688	. zdt_split, !! ptimestep
689	. zplev_omp, !! pplev
690	. zplay_omp, !! pplay
691	. zphi_omp, !! pphi
692	. zufi_omp, !! pu
693	. zvfi_omp, !! pv
694	. ztfi_omp, !! pt
695	. zqfi_omp, !! pq
696	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
697	. zdufi_omp, !! pdu
698	. zdvfi_omp, !! pdv
699	. zdtfi_omp, !! pdt
700	. zdqfi_omp, !! pdq
701	. zdpsrf_omp, !! pdpsrf
702	. tracerdyn) !! tracerdyn <-- utilite ???
703
704	endif ! of if (planet_type=="earth")
705
706
707	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
708	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
709	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
710	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
711
712	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
713	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
714	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
715	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
716
717	enddo
718
719	#endif
720	! of #ifdef CPP_PHYS
721
722	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
723	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
724	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
725	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
726	c$OMP BARRIER
727
728	do l=1,llm+1
729	do i=1,klon
730	zplev(offset+i,l)=zplev_omp(i,l)
731	enddo
732	enddo
733
734	do l=1,llm
735	do i=1,klon
736	zplay(offset+i,l)=zplay_omp(i,l)
737	enddo
738	enddo
739
740	do l=1,llm
741	do i=1,klon
742	zphi(offset+i,l)=zphi_omp(i,l)
743	enddo
744	enddo
745
746
747	do i=1,klon
748	zphis(offset+i)=zphis_omp(i)
749	enddo
750
751
752	do l=1,llm
753	presnivs(l)=presnivs_omp(l)
754	enddo
755
756	do l=1,llm
757	do i=1,klon
758	zufi(offset+i,l)=zufi_omp(i,l)
759	enddo
760	enddo
761
762	do l=1,llm
763	do i=1,klon
764	zvfi(offset+i,l)=zvfi_omp(i,l)
765	enddo
766	enddo
767
768	do l=1,llm
769	do i=1,klon
770	ztfi(offset+i,l)=ztfi_omp(i,l)
771	enddo
772	enddo
773
774	do iq=1,nqtot
775	do l=1,llm
776	do i=1,klon
777	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
778	enddo
779	enddo
780	enddo
781
782	do l=1,llm
783	do i=1,klon
784	zdufi(offset+i,l)=zdufi_omp(i,l)
785	enddo
786	enddo
787
788	do l=1,llm
789	do i=1,klon
790	zdvfi(offset+i,l)=zdvfi_omp(i,l)
791	enddo
792	enddo
793
794	do l=1,llm
795	do i=1,klon
796	zdtfi(offset+i,l)=zdtfi_omp(i,l)
797	enddo
798	enddo
799
800	do iq=1,nqtot
801	do l=1,llm
802	do i=1,klon
803	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
804	enddo
805	enddo
806	enddo
807
808	do i=1,klon
809	zdpsrf(offset+i)=zdpsrf_omp(i)
810	enddo
811
812
813	klon=klon_mpi
814	500 CONTINUE
815	c$OMP BARRIER
816
817	c$OMP MASTER
818	call stop_timer(timer_physic)
819	c$OMP END MASTER
820
821	IF (using_mpi) THEN
822
823	if (MPI_rank>0) then
824
825	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
826	DO l=1,llm
827	du_send(1:iim,l)=zdufi(1:iim,l)
828	dv_send(1:iim,l)=zdvfi(1:iim,l)
829	ENDDO
830	c$OMP END DO NOWAIT
831
832	c$OMP BARRIER
833	#ifdef CPP_MPI
834	c$OMP MASTER
835	!$OMP CRITICAL (MPI)
836	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
837	& COMM_LMDZ,Req(1),ierr)
838	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
839	& COMM_LMDZ,Req(2),ierr)
840	!$OMP END CRITICAL (MPI)
841	c$OMP END MASTER
842	#endif
843	c$OMP BARRIER
844
845	endif
846
847	if (MPI_rank<MPI_Size-1) then
848	c$OMP BARRIER
849	#ifdef CPP_MPI
850	c$OMP MASTER
851	!$OMP CRITICAL (MPI)
852	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
853	& COMM_LMDZ,Req(3),ierr)
854	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
855	& COMM_LMDZ,Req(4),ierr)
856	!$OMP END CRITICAL (MPI)
857	c$OMP END MASTER
858	#endif
859	endif
860
861	c$OMP BARRIER
862
863
864	#ifdef CPP_MPI
865	c$OMP MASTER
866	!$OMP CRITICAL (MPI)
867	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
868	call MPI_WAITALL(4,Req(1),Status,ierr)
869	else if (MPI_rank>0) then
870	call MPI_WAITALL(2,Req(1),Status,ierr)
871	else if (MPI_rank <MPI_Size-1) then
872	call MPI_WAITALL(2,Req(3),Status,ierr)
873	endif
874	!$OMP END CRITICAL (MPI)
875	c$OMP END MASTER
876	#endif
877
878	c$OMP BARRIER
879
880	ENDIF ! using_mpi
881
882
883	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
884	DO l=1,llm
885
886	zdufi2(1:klon,l)=zdufi(1:klon,l)
887	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
888
889	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
890	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
891
892	pdhfi(:,jj_begin,l)=0
893	pdqfi(:,jj_begin,l,:)=0
894	pdufi(:,jj_begin,l)=0
895	pdvfi(:,jj_begin,l)=0
896
897	if (.not. is_south_pole) then
898	pdhfi(:,jj_end,l)=0
899	pdqfi(:,jj_end,l,:)=0
900	pdufi(:,jj_end,l)=0
901	pdvfi(:,jj_end,l)=0
902	endif
903
904	ENDDO
905	c$OMP END DO NOWAIT
906
907	c$OMP MASTER
908	pdpsfi(:,jj_begin)=0
909	if (.not. is_south_pole) then
910	pdpsfi(:,jj_end)=0
911	endif
912	c$OMP END MASTER
913	c-----------------------------------------------------------------------
914	c transformation des tendances physiques en tendances dynamiques:
915	c ---------------------------------------------------------------
916
917	c tendance sur la pression :
918	c -----------------------------------
919	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
920	c
921	c 62. enthalpie potentielle
922	c ---------------------
923
924	kstart=1
925	kend=klon
926
927	if (is_north_pole) kstart=2
928	if (is_south_pole) kend=klon-1
929
930	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
931	DO l=1,llm
932
933	!CDIR ON_ADB(index_i)
934	!CDIR ON_ADB(index_j)
935	!cdir NODEP
936	do ig0=kstart,kend
937	i=index_i(ig0)
938	j=index_j(ig0)
939	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
940	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
941	enddo
942
943	if (is_north_pole) then
944	DO i=1,iip1
945	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
946	enddo
947	endif
948
949	if (is_south_pole) then
950	DO i=1,iip1
951	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
952	ENDDO
953	endif
954	ENDDO
955	c$OMP END DO NOWAIT
956
957	c 62. humidite specifique
958	c ---------------------
959	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
960	! DO iq=1,nqtot
961	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
962	! DO l=1,llm
963	!!!cdir NODEP
964	! do ig0=kstart,kend
965	! i=index_i(ig0)
966	! j=index_j(ig0)
967	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
968	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
969	! enddo
970	!
971	! if (is_north_pole) then
972	! do i=1,iip1
973	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
974	! enddo
975	! endif
976	!
977	! if (is_south_pole) then
978	! do i=1,iip1
979	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
980	! enddo
981	! endif
982	! ENDDO
983	!c$OMP END DO NOWAIT
984	! ENDDO
985
986	c 63. traceurs
987	c ------------
988	C initialisation des tendances
989
990	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
991	DO l=1,llm
992	pdqfi(:,:,l,:)=0.
993	ENDDO
994	c$OMP END DO NOWAIT
995
996	C
997	!cdir NODEP
998	DO iq=1,nqtot
999	iiq=niadv(iq)
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	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1009	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1010	ENDDO
1011
1012	IF (is_north_pole) then
1013	DO i=1,iip1
1014	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
1015	ENDDO
1016	ENDIF
1017
1018	IF (is_south_pole) then
1019	DO i=1,iip1
1020	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
1021	ENDDO
1022	ENDIF
1023
1024	ENDDO
1025	c$OMP END DO NOWAIT
1026	ENDDO
1027
1028	c 65. champ u:
1029	c ------------
1030	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1031	DO l=1,llm
1032	!CDIR ON_ADB(index_i)
1033	!CDIR ON_ADB(index_j)
1034	!cdir NODEP
1035	do ig0=kstart,kend
1036	i=index_i(ig0)
1037	j=index_j(ig0)
1038
1039	if (i/=iim) then
1040	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1041	endif
1042
1043	if (i==1) then
1044	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1045	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1046	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1047	endif
1048
1049	enddo
1050
1051	if (is_north_pole) then
1052	DO i=1,iip1
1053	pdufi(i,1,l) = 0.
1054	ENDDO
1055	endif
1056
1057	if (is_south_pole) then
1058	DO i=1,iip1
1059	pdufi(i,jjp1,l) = 0.
1060	ENDDO
1061	endif
1062
1063	ENDDO
1064	c$OMP END DO NOWAIT
1065
1066	c 67. champ v:
1067	c ------------
1068
1069	kstart=1
1070	kend=klon
1071
1072	if (is_north_pole) kstart=2
1073	if (is_south_pole) kend=klon-1-iim
1074
1075	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1076	DO l=1,llm
1077	!CDIR ON_ADB(index_i)
1078	!CDIR ON_ADB(index_j)
1079	!cdir NODEP
1080	do ig0=kstart,kend
1081	i=index_i(ig0)
1082	j=index_j(ig0)
1083	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1084	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1085	$ zdvfi2(ig0+iim,l))
1086	$ *cv(i,j)
1087	enddo
1088
1089	ENDDO
1090	c$OMP END DO NOWAIT
1091
1092
1093	c 68. champ v pres des poles:
1094	c ---------------------------
1095	c v = U * cos(long) + V * SIN(long)
1096
1097	if (is_north_pole) then
1098
1099	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1100	DO l=1,llm
1101
1102	DO i=1,iim
1103	pdvfi(i,1,l)=
1104	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1105
1106	pdvfi(i,1,l)=
1107	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1108	ENDDO
1109
1110	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1111
1112	ENDDO
1113	c$OMP END DO NOWAIT
1114
1115	endif
1116
1117	if (is_south_pole) then
1118
1119	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1120	DO l=1,llm
1121
1122	DO i=1,iim
1123	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1124	$ +zdvfi(klon,l)*SIN(rlonv(i))
1125
1126	pdvfi(i,jjm,l)=
1127	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1128	ENDDO
1129
1130	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1131
1132	ENDDO
1133	c$OMP END DO NOWAIT
1134
1135	endif
1136	c-----------------------------------------------------------------------
1137
1138	700 CONTINUE
1139
1140	firstcal = .FALSE.
1141
1142	#else
1143	write(lunout,*)
1144	& "calfis_p: for now can only work with parallel physics"
1145	stop
1146	#endif
1147	! of #ifdef CPP_PHYS
1148	RETURN
1149	END

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

Download in other formats:

Original Format