Context Navigation

calfis_p.F @ 2024

Last change on this file since 2024 was 1999, checked in by Laurent Fairhead, 12 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 30.3 KB

Line
1	!
2	! $Id: calfis_p.F 1999 2014-03-20 09:57:19Z fhourdin $
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_lmdz, ONLY : omp_chunk, using_mpi
37	USE Write_Field
38	Use Write_field_p
39	USE Times
40	USE infotrac, ONLY: nqtot, niadv, tname
41	USE control_mod, ONLY: planet_type, nsplit_phys
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,INTENT(IN) :: lafin ! .true. for the very last call to physics
115	REAL,INTENT(IN) :: jD_cur, jH_cur
116	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
117	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
118	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
119	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
120	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
121	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
122	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
123
124	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov ! not used
125	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
126	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
127	! NB: pdteta is used only to compute pcvgt which is in fact not used...
128	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
129	! NB: pdq is only used to compute pcvgq which is in fact not used...
130
131	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
132	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
133	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
134	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on dynamics grid
135
136	! tendencies (in */s) from the physics
137	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
138	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
139	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
140	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
141	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
142
143	INTEGER,PARAMETER :: longcles = 20
144	REAL,INTENT(IN) :: clesphy0( longcles ) ! unused
145
146	#ifdef CPP_PHYS
147	! Ehouarn: for now calfis_p needs some informations from physics to compile
148	c Local variables :
149	c -----------------
150
151	INTEGER i,j,l,ig0,ig,iq,iiq
152	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
153	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
154	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
155	c
156	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
157	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
158	c
159	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
160	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
161	c
162	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
163	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
164	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
165	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
166
167	c
168	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
169	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
170	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
171	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
172	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
173	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
174	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
175	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
176	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
177	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
178	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
180	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
181	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
182	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
183
184	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185	! Introduction du splitting (FH)
186	! Question pour Yann :
187	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
188	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
189	! soit allocatable (plutot par exemple que de passer une dimension
190	! dépendant du process en argument des routines) et que, du coup,
191	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
192	! Tu confirmes ?
193	! J'ai suivi le même principe pour les zdufic_omp
194	! Mais c'est surement bien que tu controles.
195	!
196
197	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
198	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
199	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
200	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
201	REAL jH_cur_split,zdt_split
202	LOGICAL debut_split,lafin_split
203	INTEGER isplit
204	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
205
206	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
207	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
208	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
209	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
210	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
211
212	LOGICAL,SAVE :: first_omp=.true.
213	c$OMP THREADPRIVATE(first_omp)
214
215	REAL zsin(iim),zcos(iim),z1(iim)
216	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
217	REAL unskap, pksurcp
218	c
219	cIM diagnostique PVteta, Amip2
220	INTEGER,PARAMETER :: ntetaSTD=3
221	REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
222	REAL PVteta(klon,ntetaSTD)
223
224	REAL SSUM
225
226	LOGICAL,SAVE :: firstcal=.true., debut=.true.
227	c$OMP THREADPRIVATE(firstcal,debut)
228
229	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
230	INTEGER :: ierr
231	#ifdef CPP_MPI
232	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
233	#else
234	INTEGER,dimension(1,4) :: Status
235	#endif
236	INTEGER, dimension(4) :: Req
237	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
238	integer :: k,kstart,kend
239	INTEGER :: offset
240
241	LOGICAL tracerdyn
242	c
243	c-----------------------------------------------------------------------
244	c
245	c 1. Initialisations :
246	c --------------------
247	c
248
249	klon=klon_mpi
250
251	PVteta(:,:)=0.
252
253	c
254	IF ( firstcal ) THEN
255	debut = .TRUE.
256	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
257	write(lunout,*) 'STOP dans calfis'
258	write(lunout,*)
259	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
260	write(lunout,*) ' ngridmx jjm iim '
261	write(lunout,*) ngridmx,jjm,iim
262	STOP
263	ENDIF
264	c$OMP MASTER
265	ALLOCATE(zpsrf(klon))
266	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
267	ALLOCATE(zphi(klon,llm),zphis(klon))
268	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
269	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
270	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
271	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
272	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
273	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
274	ALLOCATE(zdpsrf(klon))
275	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
276	ALLOCATE(flxwfi(klon,llm))
277	c$OMP END MASTER
278	c$OMP BARRIER
279	ELSE
280	debut = .FALSE.
281	ENDIF
282
283	c
284	c
285	c-----------------------------------------------------------------------
286	c 40. transformation des variables dynamiques en variables physiques:
287	c ---------------------------------------------------------------
288
289	c 41. pressions au sol (en Pascals)
290	c ----------------------------------
291
292	c$OMP MASTER
293	call start_timer(timer_physic)
294	c$OMP END MASTER
295
296	c$OMP MASTER
297	!CDIR ON_ADB(index_i)
298	!CDIR ON_ADB(index_j)
299	do ig0=1,klon
300	i=index_i(ig0)
301	j=index_j(ig0)
302	zpsrf(ig0)=pps(i,j)
303	enddo
304	c$OMP END MASTER
305
306
307	c 42. pression intercouches :
308	c
309	c -----------------------------------------------------------------
310	c .... zplev definis aux (llm +1) interfaces des couches ....
311	c .... zplay definis aux ( llm ) milieux des couches ....
312	c -----------------------------------------------------------------
313
314	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
315	c
316	unskap = 1./ kappa
317	c
318	c print *,omp_rank,'klon--->',klon
319	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
320	DO l = 1, llmp1
321	!CDIR ON_ADB(index_i)
322	!CDIR ON_ADB(index_j)
323	do ig0=1,klon
324	i=index_i(ig0)
325	j=index_j(ig0)
326	zplev( ig0,l ) = pp(i,j,l)
327	enddo
328	ENDDO
329	c$OMP END DO NOWAIT
330	c
331	c
332
333	c 43. temperature naturelle (en K) et pressions milieux couches .
334	c ---------------------------------------------------------------
335	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
336	DO l=1,llm
337	!CDIR ON_ADB(index_i)
338	!CDIR ON_ADB(index_j)
339	do ig0=1,klon
340	i=index_i(ig0)
341	j=index_j(ig0)
342	pksurcp = ppk(i,j,l) / cpp
343	zplay(ig0,l) = preff * pksurcp ** unskap
344	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
345	enddo
346
347	ENDDO
348	c$OMP END DO NOWAIT
349
350	c 43.bis traceurs
351	c ---------------
352	c
353
354	DO iq=1,nqtot
355	iiq=niadv(iq)
356	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
357	DO l=1,llm
358	!CDIR ON_ADB(index_i)
359	!CDIR ON_ADB(index_j)
360	do ig0=1,klon
361	i=index_i(ig0)
362	j=index_j(ig0)
363	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
364	enddo
365	ENDDO
366	c$OMP END DO NOWAIT
367	ENDDO
368
369
370	c Geopotentiel calcule par rapport a la surface locale:
371	c -----------------------------------------------------
372
373	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
374
375	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
376
377	c$OMP BARRIER
378
379	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
380	DO l=1,llm
381	DO ig=1,klon
382	zphi(ig,l)=zphi(ig,l)-zphis(ig)
383	ENDDO
384	ENDDO
385	c$OMP END DO NOWAIT
386
387
388	c
389	c 45. champ u:
390	c ------------
391
392	kstart=1
393	kend=klon
394
395	if (is_north_pole) kstart=2
396	if (is_south_pole) kend=klon-1
397
398	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
399	DO l=1,llm
400	!CDIR ON_ADB(index_i)
401	!CDIR ON_ADB(index_j)
402	!CDIR SPARSE
403	do ig0=kstart,kend
404	i=index_i(ig0)
405	j=index_j(ig0)
406	if (i==1) then
407	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
408	$ + pucov(1,j,l)/cu(1,j) )
409	else
410	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
411	$ + pucov(i,j,l)/cu(i,j) )
412	endif
413	enddo
414	ENDDO
415	c$OMP END DO NOWAIT
416
417	c 46.champ v:
418	c -----------
419
420	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
421	DO l=1,llm
422	!CDIR ON_ADB(index_i)
423	!CDIR ON_ADB(index_j)
424	DO ig0=kstart,kend
425	i=index_i(ig0)
426	j=index_j(ig0)
427	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
428	$ + pvcov(i,j,l)/cv(i,j) )
429
430	ENDDO
431	ENDDO
432	c$OMP END DO NOWAIT
433
434	c 47. champs de vents aux pole nord
435	c ------------------------------
436	c U = 1 / pi * integrale [ v * cos(long) * d long ]
437	c V = 1 / pi * integrale [ v * sin(long) * d long ]
438
439	if (is_north_pole) then
440	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
441	DO l=1,llm
442
443	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
444	DO i=2,iim
445	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
446	ENDDO
447
448	DO i=1,iim
449	zcos(i) = COS(rlonv(i))*z1(i)
450	zsin(i) = SIN(rlonv(i))*z1(i)
451	ENDDO
452
453	zufi(1,l) = SSUM(iim,zcos,1)/pi
454	zvfi(1,l) = SSUM(iim,zsin,1)/pi
455
456	ENDDO
457	c$OMP END DO NOWAIT
458	endif
459
460
461	c 48. champs de vents aux pole sud:
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_south_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,jjm,l)/cv(1,jjm)
471	DO i=2,iim
472	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
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(klon,l) = SSUM(iim,zcos,1)/pi
481	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
482	ENDDO
483	c$OMP END DO NOWAIT
484	endif
485
486
487	IF (is_sequential.and.(planet_type=="earth")) THEN
488	#ifdef CPP_PHYS
489	! PVtheta calls tetalevel, which is in the physics
490	cIM calcul PV a teta=350, 380, 405K
491	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
492	$ ztfi,zplay,zplev,
493	$ ntetaSTD,rtetaSTD,PVteta)
494	c
495	#endif
496	ENDIF
497
498	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
499	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
500
501	c-----------------------------------------------------------------------
502	c Appel de la physique:
503	c ---------------------
504
505
506	c$OMP BARRIER
507	if (first_omp) then
508	klon=klon_omp
509
510	allocate(zplev_omp(klon,llm+1))
511	allocate(zplay_omp(klon,llm))
512	allocate(zphi_omp(klon,llm))
513	allocate(zphis_omp(klon))
514	allocate(presnivs_omp(llm))
515	allocate(zufi_omp(klon,llm))
516	allocate(zvfi_omp(klon,llm))
517	allocate(ztfi_omp(klon,llm))
518	allocate(zqfi_omp(klon,llm,nqtot))
519	allocate(zdufi_omp(klon,llm))
520	allocate(zdvfi_omp(klon,llm))
521	allocate(zdtfi_omp(klon,llm))
522	allocate(zdqfi_omp(klon,llm,nqtot))
523	allocate(zdufic_omp(klon,llm))
524	allocate(zdvfic_omp(klon,llm))
525	allocate(zdtfic_omp(klon,llm))
526	allocate(zdqfic_omp(klon,llm,nqtot))
527	allocate(zdpsrf_omp(klon))
528	allocate(flxwfi_omp(klon,llm))
529	first_omp=.false.
530	endif
531
532
533	klon=klon_omp
534	offset=klon_omp_begin-1
535
536	do l=1,llm+1
537	do i=1,klon
538	zplev_omp(i,l)=zplev(offset+i,l)
539	enddo
540	enddo
541
542	do l=1,llm
543	do i=1,klon
544	zplay_omp(i,l)=zplay(offset+i,l)
545	enddo
546	enddo
547
548	do l=1,llm
549	do i=1,klon
550	zphi_omp(i,l)=zphi(offset+i,l)
551	enddo
552	enddo
553
554	do i=1,klon
555	zphis_omp(i)=zphis(offset+i)
556	enddo
557
558
559	do l=1,llm
560	presnivs_omp(l)=presnivs(l)
561	enddo
562
563	do l=1,llm
564	do i=1,klon
565	zufi_omp(i,l)=zufi(offset+i,l)
566	enddo
567	enddo
568
569	do l=1,llm
570	do i=1,klon
571	zvfi_omp(i,l)=zvfi(offset+i,l)
572	enddo
573	enddo
574
575	do l=1,llm
576	do i=1,klon
577	ztfi_omp(i,l)=ztfi(offset+i,l)
578	enddo
579	enddo
580
581	do iq=1,nqtot
582	do l=1,llm
583	do i=1,klon
584	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
585	enddo
586	enddo
587	enddo
588
589	do l=1,llm
590	do i=1,klon
591	zdufi_omp(i,l)=zdufi(offset+i,l)
592	enddo
593	enddo
594
595	do l=1,llm
596	do i=1,klon
597	zdvfi_omp(i,l)=zdvfi(offset+i,l)
598	enddo
599	enddo
600
601	do l=1,llm
602	do i=1,klon
603	zdtfi_omp(i,l)=zdtfi(offset+i,l)
604	enddo
605	enddo
606
607	do iq=1,nqtot
608	do l=1,llm
609	do i=1,klon
610	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
611	enddo
612	enddo
613	enddo
614
615	do i=1,klon
616	zdpsrf_omp(i)=zdpsrf(offset+i)
617	enddo
618
619	do l=1,llm
620	do i=1,klon
621	flxwfi_omp(i,l)=flxwfi(offset+i,l)
622	enddo
623	enddo
624
625	c$OMP BARRIER
626
627	!$OMP MASTER
628	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
629	!$OMP END MASTER
630	zdt_split=dtphys/nsplit_phys
631	zdufic_omp(:,:)=0.
632	zdvfic_omp(:,:)=0.
633	zdtfic_omp(:,:)=0.
634	zdqfic_omp(:,:,:)=0.
635
636	#ifdef CPP_PHYS
637	do isplit=1,nsplit_phys
638
639	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
640	debut_split=debut.and.isplit==1
641	lafin_split=lafin.and.isplit==nsplit_phys
642
643	if (planet_type=="earth") then
644
645	CALL physiq (klon,
646	. llm,
647	. debut_split,
648	. lafin_split,
649	. jD_cur,
650	. jH_cur_split,
651	. zdt_split,
652	. zplev_omp,
653	. zplay_omp,
654	. zphi_omp,
655	. zphis_omp,
656	. presnivs_omp,
657	. clesphy0,
658	. zufi_omp,
659	. zvfi_omp,
660	. ztfi_omp,
661	. zqfi_omp,
662	c#ifdef INCA
663	. flxwfi_omp,
664	c#endif
665	. zdufi_omp,
666	. zdvfi_omp,
667	. zdtfi_omp,
668	. zdqfi_omp,
669	. zdpsrf_omp,
670	cIM diagnostique PVteta, Amip2
671	. pducov,
672	. PVteta)
673
674	else if ( planet_type=="generic" ) then
675
676	CALL physiq (klon, !! ngrid
677	. llm, !! nlayer
678	. nqtot, !! nq
679	. tname, !! tracer names from dynamical core (given in infotrac)
680	. debut_split, !! firstcall
681	. lafin_split, !! lastcall
682	. jD_cur, !! pday. see leapfrog_p
683	. jH_cur_split, !! ptime "fraction of day"
684	. zdt_split, !! ptimestep
685	. zplev_omp, !! pplev
686	. zplay_omp, !! pplay
687	. zphi_omp, !! pphi
688	. zufi_omp, !! pu
689	. zvfi_omp, !! pv
690	. ztfi_omp, !! pt
691	. zqfi_omp, !! pq
692	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
693	. zdufi_omp, !! pdu
694	. zdvfi_omp, !! pdv
695	. zdtfi_omp, !! pdt
696	. zdqfi_omp, !! pdq
697	. zdpsrf_omp, !! pdpsrf
698	. tracerdyn) !! tracerdyn <-- utilite ???
699
700	endif ! of if (planet_type=="earth")
701
702
703	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
704	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
705	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
706	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
707
708	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
709	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
710	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
711	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
712
713	enddo
714
715	#endif
716	! of #ifdef CPP_PHYS
717
718	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
719	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
720	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
721	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
722	c$OMP BARRIER
723
724	do l=1,llm+1
725	do i=1,klon
726	zplev(offset+i,l)=zplev_omp(i,l)
727	enddo
728	enddo
729
730	do l=1,llm
731	do i=1,klon
732	zplay(offset+i,l)=zplay_omp(i,l)
733	enddo
734	enddo
735
736	do l=1,llm
737	do i=1,klon
738	zphi(offset+i,l)=zphi_omp(i,l)
739	enddo
740	enddo
741
742
743	do i=1,klon
744	zphis(offset+i)=zphis_omp(i)
745	enddo
746
747
748	do l=1,llm
749	presnivs(l)=presnivs_omp(l)
750	enddo
751
752	do l=1,llm
753	do i=1,klon
754	zufi(offset+i,l)=zufi_omp(i,l)
755	enddo
756	enddo
757
758	do l=1,llm
759	do i=1,klon
760	zvfi(offset+i,l)=zvfi_omp(i,l)
761	enddo
762	enddo
763
764	do l=1,llm
765	do i=1,klon
766	ztfi(offset+i,l)=ztfi_omp(i,l)
767	enddo
768	enddo
769
770	do iq=1,nqtot
771	do l=1,llm
772	do i=1,klon
773	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
774	enddo
775	enddo
776	enddo
777
778	do l=1,llm
779	do i=1,klon
780	zdufi(offset+i,l)=zdufi_omp(i,l)
781	enddo
782	enddo
783
784	do l=1,llm
785	do i=1,klon
786	zdvfi(offset+i,l)=zdvfi_omp(i,l)
787	enddo
788	enddo
789
790	do l=1,llm
791	do i=1,klon
792	zdtfi(offset+i,l)=zdtfi_omp(i,l)
793	enddo
794	enddo
795
796	do iq=1,nqtot
797	do l=1,llm
798	do i=1,klon
799	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
800	enddo
801	enddo
802	enddo
803
804	do i=1,klon
805	zdpsrf(offset+i)=zdpsrf_omp(i)
806	enddo
807
808
809	klon=klon_mpi
810	500 CONTINUE
811	c$OMP BARRIER
812
813	c$OMP MASTER
814	call stop_timer(timer_physic)
815	c$OMP END MASTER
816
817	IF (using_mpi) THEN
818
819	if (MPI_rank>0) then
820
821	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
822	DO l=1,llm
823	du_send(1:iim,l)=zdufi(1:iim,l)
824	dv_send(1:iim,l)=zdvfi(1:iim,l)
825	ENDDO
826	c$OMP END DO NOWAIT
827
828	c$OMP BARRIER
829	#ifdef CPP_MPI
830	c$OMP MASTER
831	!$OMP CRITICAL (MPI)
832	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
833	& COMM_LMDZ,Req(1),ierr)
834	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
835	& COMM_LMDZ,Req(2),ierr)
836	!$OMP END CRITICAL (MPI)
837	c$OMP END MASTER
838	#endif
839	c$OMP BARRIER
840
841	endif
842
843	if (MPI_rank<MPI_Size-1) then
844	c$OMP BARRIER
845	#ifdef CPP_MPI
846	c$OMP MASTER
847	!$OMP CRITICAL (MPI)
848	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
849	& COMM_LMDZ,Req(3),ierr)
850	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
851	& COMM_LMDZ,Req(4),ierr)
852	!$OMP END CRITICAL (MPI)
853	c$OMP END MASTER
854	#endif
855	endif
856
857	c$OMP BARRIER
858
859
860	#ifdef CPP_MPI
861	c$OMP MASTER
862	!$OMP CRITICAL (MPI)
863	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
864	call MPI_WAITALL(4,Req(1),Status,ierr)
865	else if (MPI_rank>0) then
866	call MPI_WAITALL(2,Req(1),Status,ierr)
867	else if (MPI_rank <MPI_Size-1) then
868	call MPI_WAITALL(2,Req(3),Status,ierr)
869	endif
870	!$OMP END CRITICAL (MPI)
871	c$OMP END MASTER
872	#endif
873
874	c$OMP BARRIER
875
876	ENDIF ! using_mpi
877
878
879	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
880	DO l=1,llm
881
882	zdufi2(1:klon,l)=zdufi(1:klon,l)
883	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
884
885	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
886	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
887
888	pdhfi(:,jj_begin,l)=0
889	pdqfi(:,jj_begin,l,:)=0
890	pdufi(:,jj_begin,l)=0
891	pdvfi(:,jj_begin,l)=0
892
893	if (.not. is_south_pole) then
894	pdhfi(:,jj_end,l)=0
895	pdqfi(:,jj_end,l,:)=0
896	pdufi(:,jj_end,l)=0
897	pdvfi(:,jj_end,l)=0
898	endif
899
900	ENDDO
901	c$OMP END DO NOWAIT
902
903	c$OMP MASTER
904	pdpsfi(:,jj_begin)=0
905	if (.not. is_south_pole) then
906	pdpsfi(:,jj_end)=0
907	endif
908	c$OMP END MASTER
909	c-----------------------------------------------------------------------
910	c transformation des tendances physiques en tendances dynamiques:
911	c ---------------------------------------------------------------
912
913	c tendance sur la pression :
914	c -----------------------------------
915	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
916	c
917	c 62. enthalpie potentielle
918	c ---------------------
919
920	kstart=1
921	kend=klon
922
923	if (is_north_pole) kstart=2
924	if (is_south_pole) kend=klon-1
925
926	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
927	DO l=1,llm
928
929	!CDIR ON_ADB(index_i)
930	!CDIR ON_ADB(index_j)
931	!cdir NODEP
932	do ig0=kstart,kend
933	i=index_i(ig0)
934	j=index_j(ig0)
935	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
936	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
937	enddo
938
939	if (is_north_pole) then
940	DO i=1,iip1
941	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
942	enddo
943	endif
944
945	if (is_south_pole) then
946	DO i=1,iip1
947	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
948	ENDDO
949	endif
950	ENDDO
951	c$OMP END DO NOWAIT
952
953	c 62. humidite specifique
954	c ---------------------
955	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
956	! DO iq=1,nqtot
957	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
958	! DO l=1,llm
959	!!!cdir NODEP
960	! do ig0=kstart,kend
961	! i=index_i(ig0)
962	! j=index_j(ig0)
963	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
964	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
965	! enddo
966	!
967	! if (is_north_pole) then
968	! do i=1,iip1
969	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
970	! enddo
971	! endif
972	!
973	! if (is_south_pole) then
974	! do i=1,iip1
975	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
976	! enddo
977	! endif
978	! ENDDO
979	!c$OMP END DO NOWAIT
980	! ENDDO
981
982	c 63. traceurs
983	c ------------
984	C initialisation des tendances
985
986	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
987	DO l=1,llm
988	pdqfi(:,:,l,:)=0.
989	ENDDO
990	c$OMP END DO NOWAIT
991
992	C
993	!cdir NODEP
994	DO iq=1,nqtot
995	iiq=niadv(iq)
996	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
997	DO l=1,llm
998	!CDIR ON_ADB(index_i)
999	!CDIR ON_ADB(index_j)
1000	!cdir NODEP
1001	DO ig0=kstart,kend
1002	i=index_i(ig0)
1003	j=index_j(ig0)
1004	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
1005	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
1006	ENDDO
1007
1008	IF (is_north_pole) then
1009	DO i=1,iip1
1010	pdqfi(i,1,l,iiq) = 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,iiq) = zdqfi(klon,l,iq)
1017	ENDDO
1018	ENDIF
1019
1020	ENDDO
1021	c$OMP END DO NOWAIT
1022	ENDDO
1023
1024	c 65. champ u:
1025	c ------------
1026	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1027	DO l=1,llm
1028	!CDIR ON_ADB(index_i)
1029	!CDIR ON_ADB(index_j)
1030	!cdir NODEP
1031	do ig0=kstart,kend
1032	i=index_i(ig0)
1033	j=index_j(ig0)
1034
1035	if (i/=iim) then
1036	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1037	endif
1038
1039	if (i==1) then
1040	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1041	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1042	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1043	endif
1044
1045	enddo
1046
1047	if (is_north_pole) then
1048	DO i=1,iip1
1049	pdufi(i,1,l) = 0.
1050	ENDDO
1051	endif
1052
1053	if (is_south_pole) then
1054	DO i=1,iip1
1055	pdufi(i,jjp1,l) = 0.
1056	ENDDO
1057	endif
1058
1059	ENDDO
1060	c$OMP END DO NOWAIT
1061
1062	c 67. champ v:
1063	c ------------
1064
1065	kstart=1
1066	kend=klon
1067
1068	if (is_north_pole) kstart=2
1069	if (is_south_pole) kend=klon-1-iim
1070
1071	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1072	DO l=1,llm
1073	!CDIR ON_ADB(index_i)
1074	!CDIR ON_ADB(index_j)
1075	!cdir NODEP
1076	do ig0=kstart,kend
1077	i=index_i(ig0)
1078	j=index_j(ig0)
1079	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1080	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1081	$ zdvfi2(ig0+iim,l))
1082	$ *cv(i,j)
1083	enddo
1084
1085	ENDDO
1086	c$OMP END DO NOWAIT
1087
1088
1089	c 68. champ v pres des poles:
1090	c ---------------------------
1091	c v = U * cos(long) + V * SIN(long)
1092
1093	if (is_north_pole) then
1094
1095	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1096	DO l=1,llm
1097
1098	DO i=1,iim
1099	pdvfi(i,1,l)=
1100	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1101
1102	pdvfi(i,1,l)=
1103	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1104	ENDDO
1105
1106	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1107
1108	ENDDO
1109	c$OMP END DO NOWAIT
1110
1111	endif
1112
1113	if (is_south_pole) then
1114
1115	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1116	DO l=1,llm
1117
1118	DO i=1,iim
1119	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1120	$ +zdvfi(klon,l)*SIN(rlonv(i))
1121
1122	pdvfi(i,jjm,l)=
1123	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1124	ENDDO
1125
1126	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1127
1128	ENDDO
1129	c$OMP END DO NOWAIT
1130
1131	endif
1132	c-----------------------------------------------------------------------
1133
1134	700 CONTINUE
1135
1136	firstcal = .FALSE.
1137
1138	#else
1139	write(lunout,*)
1140	& "calfis_p: for now can only work with parallel physics"
1141	stop
1142	#endif
1143	! of #ifdef CPP_PHYS
1144	RETURN
1145	END

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

Download in other formats:

Original Format