Context Navigation

calfis_p.F @ 2239

Last change on this file since 2239 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

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: 29.6 KB

Line
1	!
2	! $Id: calfis_p.F 2239 2015-03-23 07:27:30Z emillour $
3	!
4	C
5	C
6	SUBROUTINE calfis_p(lafin,
7	$ jD_cur, jH_cur,
8	$ pucov,
9	$ pvcov,
10	$ pteta,
11	$ pq,
12	$ pmasse,
13	$ pps,
14	$ pp,
15	$ ppk,
16	$ pphis,
17	$ pphi,
18	$ pducov,
19	$ pdvcov,
20	$ pdteta,
21	$ pdq,
22	$ flxw,
23	$ pdufi,
24	$ pdvfi,
25	$ pdhfi,
26	$ pdqfi,
27	$ pdpsfi)
28	#ifdef CPP_PHYS
29	! If using physics
30	USE dimphy
31	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
32	USE mod_interface_dyn_phys
33	USE IOPHY
34	#endif
35	#ifdef CPP_PARA
36	USE parallel_lmdz, ONLY : omp_chunk, using_mpi
37	USE Write_Field
38	Use Write_field_p
39	USE Times
40	#endif
41	USE infotrac, ONLY: nqtot, niadv, tname
42	USE control_mod, ONLY: planet_type, nsplit_phys
43
44	IMPLICIT NONE
45	c=======================================================================
46	c
47	c 1. rearrangement des tableaux et transformation
48	c variables dynamiques > variables physiques
49	c 2. calcul des termes physiques
50	c 3. retransformation des tendances physiques en tendances dynamiques
51	c
52	c remarques:
53	c ----------
54	c
55	c - les vents sont donnes dans la physique par leurs composantes
56	c naturelles.
57	c - la variable thermodynamique de la physique est une variable
58	c intensive : T
59	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
60	c - les deux seules variables dependant de la geometrie necessaires
61	c pour la physique sont la latitude pour le rayonnement et
62	c l'aire de la maille quand on veut integrer une grandeur
63	c horizontalement.
64	c - les points de la physique sont les points scalaires de la
65	c la dynamique; numerotation:
66	c 1 pour le pole nord
67	c (jjm-1)*iim pour l'interieur du domaine
68	c ngridmx pour le pole sud
69	c ---> ngridmx=2+(jjm-1)*iim
70	c
71	c Input :
72	c -------
73	c ecritphy frequence d'ecriture (en jours)de histphy
74	c pucov covariant zonal velocity
75	c pvcov covariant meridional velocity
76	c pteta potential temperature
77	c pps surface pressure
78	c pmasse masse d'air dans chaque maille
79	c pts surface temperature (K)
80	c callrad clef d'appel au rayonnement
81	c
82	c Output :
83	c --------
84	c pdufi tendency for the natural zonal velocity (ms-1)
85	c pdvfi tendency for the natural meridional velocity
86	c pdhfi tendency for the potential temperature
87	c pdtsfi tendency for the surface temperature
88	c
89	c pdtrad radiative tendencies \ both input
90	c pfluxrad radiative fluxes / and output
91	c
92	c=======================================================================
93	c
94	c-----------------------------------------------------------------------
95	c
96	c 0. Declarations :
97	c ------------------
98
99	#include "dimensions.h"
100	#include "paramet.h"
101	#include "temps.h"
102
103	INTEGER ngridmx
104	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
105
106	#include "comconst.h"
107	#include "comvert.h"
108	#include "comgeom2.h"
109	#include "iniprint.h"
110	#ifdef CPP_MPI
111	include 'mpif.h'
112	#endif
113	c Arguments :
114	c -----------
115	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
116	REAL,INTENT(IN) :: jD_cur, jH_cur
117	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
118	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
119	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
120	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
121	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
122	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
123	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
124
125	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov ! not used
126	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
127	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
128	! NB: pdteta is used only to compute pcvgt which is in fact not used...
129	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
130	! NB: pdq is only used to compute pcvgq which is in fact not used...
131
132	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
133	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
134	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
135	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on dynamics grid
136
137	! tendencies (in */s) from the physics
138	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
139	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
140	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
141	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
142	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
143
144	#ifdef CPP_PARA
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	REAL SSUM
219
220	LOGICAL,SAVE :: firstcal=.true., debut=.true.
221	c$OMP THREADPRIVATE(firstcal,debut)
222
223	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
224	INTEGER :: ierr
225	#ifdef CPP_MPI
226	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
227	#else
228	INTEGER,dimension(1,4) :: Status
229	#endif
230	INTEGER, dimension(4) :: Req
231	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
232	integer :: k,kstart,kend
233	INTEGER :: offset
234
235	LOGICAL tracerdyn
236	c
237	c-----------------------------------------------------------------------
238	c
239	c 1. Initialisations :
240	c --------------------
241	c
242
243	klon=klon_mpi
244
245	c
246	IF ( firstcal ) THEN
247	debut = .TRUE.
248	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
249	write(lunout,*) 'STOP dans calfis'
250	write(lunout,*)
251	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
252	write(lunout,*) ' ngridmx jjm iim '
253	write(lunout,*) ngridmx,jjm,iim
254	STOP
255	ENDIF
256	c$OMP MASTER
257	ALLOCATE(zpsrf(klon))
258	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
259	ALLOCATE(zphi(klon,llm),zphis(klon))
260	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
261	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
262	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
263	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
264	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
265	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
266	ALLOCATE(zdpsrf(klon))
267	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
268	ALLOCATE(flxwfi(klon,llm))
269	c$OMP END MASTER
270	c$OMP BARRIER
271	ELSE
272	debut = .FALSE.
273	ENDIF
274
275	c
276	c
277	c-----------------------------------------------------------------------
278	c 40. transformation des variables dynamiques en variables physiques:
279	c ---------------------------------------------------------------
280
281	c 41. pressions au sol (en Pascals)
282	c ----------------------------------
283
284	c$OMP MASTER
285	call start_timer(timer_physic)
286	c$OMP END MASTER
287
288	c$OMP MASTER
289	!CDIR ON_ADB(index_i)
290	!CDIR ON_ADB(index_j)
291	do ig0=1,klon
292	i=index_i(ig0)
293	j=index_j(ig0)
294	zpsrf(ig0)=pps(i,j)
295	enddo
296	c$OMP END MASTER
297
298
299	c 42. pression intercouches :
300	c
301	c -----------------------------------------------------------------
302	c .... zplev definis aux (llm +1) interfaces des couches ....
303	c .... zplay definis aux ( llm ) milieux des couches ....
304	c -----------------------------------------------------------------
305
306	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
307	c
308	unskap = 1./ kappa
309	c
310	c print *,omp_rank,'klon--->',klon
311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312	DO l = 1, llmp1
313	!CDIR ON_ADB(index_i)
314	!CDIR ON_ADB(index_j)
315	do ig0=1,klon
316	i=index_i(ig0)
317	j=index_j(ig0)
318	zplev( ig0,l ) = pp(i,j,l)
319	enddo
320	ENDDO
321	c$OMP END DO NOWAIT
322	c
323	c
324
325	c 43. temperature naturelle (en K) et pressions milieux couches .
326	c ---------------------------------------------------------------
327	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
328	DO l=1,llm
329	!CDIR ON_ADB(index_i)
330	!CDIR ON_ADB(index_j)
331	do ig0=1,klon
332	i=index_i(ig0)
333	j=index_j(ig0)
334	pksurcp = ppk(i,j,l) / cpp
335	zplay(ig0,l) = preff * pksurcp ** unskap
336	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
337	enddo
338
339	ENDDO
340	c$OMP END DO NOWAIT
341
342	c 43.bis traceurs
343	c ---------------
344	c
345
346	DO iq=1,nqtot
347	iiq=niadv(iq)
348	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
349	DO l=1,llm
350	!CDIR ON_ADB(index_i)
351	!CDIR ON_ADB(index_j)
352	do ig0=1,klon
353	i=index_i(ig0)
354	j=index_j(ig0)
355	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
356	enddo
357	ENDDO
358	c$OMP END DO NOWAIT
359	ENDDO
360
361
362	c Geopotentiel calcule par rapport a la surface locale:
363	c -----------------------------------------------------
364
365	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
366
367	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
368
369	c$OMP BARRIER
370
371	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
372	DO l=1,llm
373	DO ig=1,klon
374	zphi(ig,l)=zphi(ig,l)-zphis(ig)
375	ENDDO
376	ENDDO
377	c$OMP END DO NOWAIT
378
379
380	c
381	c 45. champ u:
382	c ------------
383
384	kstart=1
385	kend=klon
386
387	if (is_north_pole) kstart=2
388	if (is_south_pole) kend=klon-1
389
390	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
391	DO l=1,llm
392	!CDIR ON_ADB(index_i)
393	!CDIR ON_ADB(index_j)
394	!CDIR SPARSE
395	do ig0=kstart,kend
396	i=index_i(ig0)
397	j=index_j(ig0)
398	if (i==1) then
399	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
400	$ + pucov(1,j,l)/cu(1,j) )
401	else
402	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
403	$ + pucov(i,j,l)/cu(i,j) )
404	endif
405	enddo
406	ENDDO
407	c$OMP END DO NOWAIT
408
409	c 46.champ v:
410	c -----------
411
412	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
413	DO l=1,llm
414	!CDIR ON_ADB(index_i)
415	!CDIR ON_ADB(index_j)
416	DO ig0=kstart,kend
417	i=index_i(ig0)
418	j=index_j(ig0)
419	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
420	$ + pvcov(i,j,l)/cv(i,j) )
421
422	ENDDO
423	ENDDO
424	c$OMP END DO NOWAIT
425
426	c 47. champs de vents aux pole nord
427	c ------------------------------
428	c U = 1 / pi * integrale [ v * cos(long) * d long ]
429	c V = 1 / pi * integrale [ v * sin(long) * d long ]
430
431	if (is_north_pole) then
432	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
433	DO l=1,llm
434
435	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
436	DO i=2,iim
437	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
438	ENDDO
439
440	DO i=1,iim
441	zcos(i) = COS(rlonv(i))*z1(i)
442	zsin(i) = SIN(rlonv(i))*z1(i)
443	ENDDO
444
445	zufi(1,l) = SSUM(iim,zcos,1)/pi
446	zvfi(1,l) = SSUM(iim,zsin,1)/pi
447
448	ENDDO
449	c$OMP END DO NOWAIT
450	endif
451
452
453	c 48. champs de vents aux pole sud:
454	c ---------------------------------
455	c U = 1 / pi * integrale [ v * cos(long) * d long ]
456	c V = 1 / pi * integrale [ v * sin(long) * d long ]
457
458	if (is_south_pole) then
459	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
460	DO l=1,llm
461
462	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
463	DO i=2,iim
464	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
465	ENDDO
466
467	DO i=1,iim
468	zcos(i) = COS(rlonv(i))*z1(i)
469	zsin(i) = SIN(rlonv(i))*z1(i)
470	ENDDO
471
472	zufi(klon,l) = SSUM(iim,zcos,1)/pi
473	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
474	ENDDO
475	c$OMP END DO NOWAIT
476	endif
477
478	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
479	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
480
481	c-----------------------------------------------------------------------
482	c Appel de la physique:
483	c ---------------------
484
485
486	c$OMP BARRIER
487	if (first_omp) then
488	klon=klon_omp
489
490	allocate(zplev_omp(klon,llm+1))
491	allocate(zplay_omp(klon,llm))
492	allocate(zphi_omp(klon,llm))
493	allocate(zphis_omp(klon))
494	allocate(presnivs_omp(llm))
495	allocate(zufi_omp(klon,llm))
496	allocate(zvfi_omp(klon,llm))
497	allocate(ztfi_omp(klon,llm))
498	allocate(zqfi_omp(klon,llm,nqtot))
499	allocate(zdufi_omp(klon,llm))
500	allocate(zdvfi_omp(klon,llm))
501	allocate(zdtfi_omp(klon,llm))
502	allocate(zdqfi_omp(klon,llm,nqtot))
503	allocate(zdufic_omp(klon,llm))
504	allocate(zdvfic_omp(klon,llm))
505	allocate(zdtfic_omp(klon,llm))
506	allocate(zdqfic_omp(klon,llm,nqtot))
507	allocate(zdpsrf_omp(klon))
508	allocate(flxwfi_omp(klon,llm))
509	first_omp=.false.
510	endif
511
512
513	klon=klon_omp
514	offset=klon_omp_begin-1
515
516	do l=1,llm+1
517	do i=1,klon
518	zplev_omp(i,l)=zplev(offset+i,l)
519	enddo
520	enddo
521
522	do l=1,llm
523	do i=1,klon
524	zplay_omp(i,l)=zplay(offset+i,l)
525	enddo
526	enddo
527
528	do l=1,llm
529	do i=1,klon
530	zphi_omp(i,l)=zphi(offset+i,l)
531	enddo
532	enddo
533
534	do i=1,klon
535	zphis_omp(i)=zphis(offset+i)
536	enddo
537
538
539	do l=1,llm
540	presnivs_omp(l)=presnivs(l)
541	enddo
542
543	do l=1,llm
544	do i=1,klon
545	zufi_omp(i,l)=zufi(offset+i,l)
546	enddo
547	enddo
548
549	do l=1,llm
550	do i=1,klon
551	zvfi_omp(i,l)=zvfi(offset+i,l)
552	enddo
553	enddo
554
555	do l=1,llm
556	do i=1,klon
557	ztfi_omp(i,l)=ztfi(offset+i,l)
558	enddo
559	enddo
560
561	do iq=1,nqtot
562	do l=1,llm
563	do i=1,klon
564	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
565	enddo
566	enddo
567	enddo
568
569	do l=1,llm
570	do i=1,klon
571	zdufi_omp(i,l)=zdufi(offset+i,l)
572	enddo
573	enddo
574
575	do l=1,llm
576	do i=1,klon
577	zdvfi_omp(i,l)=zdvfi(offset+i,l)
578	enddo
579	enddo
580
581	do l=1,llm
582	do i=1,klon
583	zdtfi_omp(i,l)=zdtfi(offset+i,l)
584	enddo
585	enddo
586
587	do iq=1,nqtot
588	do l=1,llm
589	do i=1,klon
590	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
591	enddo
592	enddo
593	enddo
594
595	do i=1,klon
596	zdpsrf_omp(i)=zdpsrf(offset+i)
597	enddo
598
599	do l=1,llm
600	do i=1,klon
601	flxwfi_omp(i,l)=flxwfi(offset+i,l)
602	enddo
603	enddo
604
605	c$OMP BARRIER
606
607	!$OMP MASTER
608	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
609	!$OMP END MASTER
610	zdt_split=dtphys/nsplit_phys
611	zdufic_omp(:,:)=0.
612	zdvfic_omp(:,:)=0.
613	zdtfic_omp(:,:)=0.
614	zdqfic_omp(:,:,:)=0.
615
616	#ifdef CPP_PHYS
617	do isplit=1,nsplit_phys
618
619	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
620	debut_split=debut.and.isplit==1
621	lafin_split=lafin.and.isplit==nsplit_phys
622
623	if (planet_type=="earth") then
624
625	CALL physiq (klon,
626	. llm,
627	. debut_split,
628	. lafin_split,
629	. jD_cur,
630	. jH_cur_split,
631	. zdt_split,
632	. zplev_omp,
633	. zplay_omp,
634	. zphi_omp,
635	. zphis_omp,
636	. presnivs_omp,
637	. zufi_omp,
638	. zvfi_omp,
639	. ztfi_omp,
640	. zqfi_omp,
641	. flxwfi_omp,
642	. zdufi_omp,
643	. zdvfi_omp,
644	. zdtfi_omp,
645	. zdqfi_omp,
646	. zdpsrf_omp,
647	. pducov)
648
649	else if ( planet_type=="generic" ) then
650
651	CALL physiq (klon, !! ngrid
652	. llm, !! nlayer
653	. nqtot, !! nq
654	. tname, !! tracer names from dynamical core (given in infotrac)
655	. debut_split, !! firstcall
656	. lafin_split, !! lastcall
657	. jD_cur, !! pday. see leapfrog_p
658	. jH_cur_split, !! ptime "fraction of day"
659	. zdt_split, !! ptimestep
660	. zplev_omp, !! pplev
661	. zplay_omp, !! pplay
662	. zphi_omp, !! pphi
663	. zufi_omp, !! pu
664	. zvfi_omp, !! pv
665	. ztfi_omp, !! pt
666	. zqfi_omp, !! pq
667	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
668	. zdufi_omp, !! pdu
669	. zdvfi_omp, !! pdv
670	. zdtfi_omp, !! pdt
671	. zdqfi_omp, !! pdq
672	. zdpsrf_omp, !! pdpsrf
673	. tracerdyn) !! tracerdyn <-- utilite ???
674
675	endif ! of if (planet_type=="earth")
676
677
678	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
679	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
680	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
681	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
682
683	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
684	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
685	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
686	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
687
688	enddo
689
690	#endif
691	! of #ifdef CPP_PHYS
692
693	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
694	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
695	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
696	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
697	c$OMP BARRIER
698
699	do l=1,llm+1
700	do i=1,klon
701	zplev(offset+i,l)=zplev_omp(i,l)
702	enddo
703	enddo
704
705	do l=1,llm
706	do i=1,klon
707	zplay(offset+i,l)=zplay_omp(i,l)
708	enddo
709	enddo
710
711	do l=1,llm
712	do i=1,klon
713	zphi(offset+i,l)=zphi_omp(i,l)
714	enddo
715	enddo
716
717
718	do i=1,klon
719	zphis(offset+i)=zphis_omp(i)
720	enddo
721
722
723	do l=1,llm
724	presnivs(l)=presnivs_omp(l)
725	enddo
726
727	do l=1,llm
728	do i=1,klon
729	zufi(offset+i,l)=zufi_omp(i,l)
730	enddo
731	enddo
732
733	do l=1,llm
734	do i=1,klon
735	zvfi(offset+i,l)=zvfi_omp(i,l)
736	enddo
737	enddo
738
739	do l=1,llm
740	do i=1,klon
741	ztfi(offset+i,l)=ztfi_omp(i,l)
742	enddo
743	enddo
744
745	do iq=1,nqtot
746	do l=1,llm
747	do i=1,klon
748	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
749	enddo
750	enddo
751	enddo
752
753	do l=1,llm
754	do i=1,klon
755	zdufi(offset+i,l)=zdufi_omp(i,l)
756	enddo
757	enddo
758
759	do l=1,llm
760	do i=1,klon
761	zdvfi(offset+i,l)=zdvfi_omp(i,l)
762	enddo
763	enddo
764
765	do l=1,llm
766	do i=1,klon
767	zdtfi(offset+i,l)=zdtfi_omp(i,l)
768	enddo
769	enddo
770
771	do iq=1,nqtot
772	do l=1,llm
773	do i=1,klon
774	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
775	enddo
776	enddo
777	enddo
778
779	do i=1,klon
780	zdpsrf(offset+i)=zdpsrf_omp(i)
781	enddo
782
783
784	klon=klon_mpi
785	500 CONTINUE
786	c$OMP BARRIER
787
788	c$OMP MASTER
789	call stop_timer(timer_physic)
790	c$OMP END MASTER
791
792	IF (using_mpi) THEN
793
794	if (MPI_rank>0) then
795
796	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
797	DO l=1,llm
798	du_send(1:iim,l)=zdufi(1:iim,l)
799	dv_send(1:iim,l)=zdvfi(1:iim,l)
800	ENDDO
801	c$OMP END DO NOWAIT
802
803	c$OMP BARRIER
804	#ifdef CPP_MPI
805	c$OMP MASTER
806	!$OMP CRITICAL (MPI)
807	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
808	& COMM_LMDZ,Req(1),ierr)
809	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
810	& COMM_LMDZ,Req(2),ierr)
811	!$OMP END CRITICAL (MPI)
812	c$OMP END MASTER
813	#endif
814	c$OMP BARRIER
815
816	endif
817
818	if (MPI_rank<MPI_Size-1) then
819	c$OMP BARRIER
820	#ifdef CPP_MPI
821	c$OMP MASTER
822	!$OMP CRITICAL (MPI)
823	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
824	& COMM_LMDZ,Req(3),ierr)
825	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
826	& COMM_LMDZ,Req(4),ierr)
827	!$OMP END CRITICAL (MPI)
828	c$OMP END MASTER
829	#endif
830	endif
831
832	c$OMP BARRIER
833
834
835	#ifdef CPP_MPI
836	c$OMP MASTER
837	!$OMP CRITICAL (MPI)
838	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
839	call MPI_WAITALL(4,Req(1),Status,ierr)
840	else if (MPI_rank>0) then
841	call MPI_WAITALL(2,Req(1),Status,ierr)
842	else if (MPI_rank <MPI_Size-1) then
843	call MPI_WAITALL(2,Req(3),Status,ierr)
844	endif
845	!$OMP END CRITICAL (MPI)
846	c$OMP END MASTER
847	#endif
848
849	c$OMP BARRIER
850
851	ENDIF ! using_mpi
852
853
854	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
855	DO l=1,llm
856
857	zdufi2(1:klon,l)=zdufi(1:klon,l)
858	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
859
860	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
861	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
862
863	pdhfi(:,jj_begin,l)=0
864	pdqfi(:,jj_begin,l,:)=0
865	pdufi(:,jj_begin,l)=0
866	pdvfi(:,jj_begin,l)=0
867
868	if (.not. is_south_pole) then
869	pdhfi(:,jj_end,l)=0
870	pdqfi(:,jj_end,l,:)=0
871	pdufi(:,jj_end,l)=0
872	pdvfi(:,jj_end,l)=0
873	endif
874
875	ENDDO
876	c$OMP END DO NOWAIT
877
878	c$OMP MASTER
879	pdpsfi(:,jj_begin)=0
880	if (.not. is_south_pole) then
881	pdpsfi(:,jj_end)=0
882	endif
883	c$OMP END MASTER
884	c-----------------------------------------------------------------------
885	c transformation des tendances physiques en tendances dynamiques:
886	c ---------------------------------------------------------------
887
888	c tendance sur la pression :
889	c -----------------------------------
890	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
891	c
892	c 62. enthalpie potentielle
893	c ---------------------
894
895	kstart=1
896	kend=klon
897
898	if (is_north_pole) kstart=2
899	if (is_south_pole) kend=klon-1
900
901	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
902	DO l=1,llm
903
904	!CDIR ON_ADB(index_i)
905	!CDIR ON_ADB(index_j)
906	!cdir NODEP
907	do ig0=kstart,kend
908	i=index_i(ig0)
909	j=index_j(ig0)
910	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
911	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
912	enddo
913
914	if (is_north_pole) then
915	DO i=1,iip1
916	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
917	enddo
918	endif
919
920	if (is_south_pole) then
921	DO i=1,iip1
922	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
923	ENDDO
924	endif
925	ENDDO
926	c$OMP END DO NOWAIT
927
928	c 62. humidite specifique
929	c ---------------------
930	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
931	! DO iq=1,nqtot
932	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
933	! DO l=1,llm
934	!!!cdir NODEP
935	! do ig0=kstart,kend
936	! i=index_i(ig0)
937	! j=index_j(ig0)
938	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
939	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
940	! enddo
941	!
942	! if (is_north_pole) then
943	! do i=1,iip1
944	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
945	! enddo
946	! endif
947	!
948	! if (is_south_pole) then
949	! do i=1,iip1
950	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
951	! enddo
952	! endif
953	! ENDDO
954	!c$OMP END DO NOWAIT
955	! ENDDO
956
957	c 63. traceurs
958	c ------------
959	C initialisation des tendances
960
961	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
962	DO l=1,llm
963	pdqfi(:,:,l,:)=0.
964	ENDDO
965	c$OMP END DO NOWAIT
966
967	C
968	!cdir NODEP
969	DO iq=1,nqtot
970	iiq=niadv(iq)
971	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
972	DO l=1,llm
973	!CDIR ON_ADB(index_i)
974	!CDIR ON_ADB(index_j)
975	!cdir NODEP
976	DO ig0=kstart,kend
977	i=index_i(ig0)
978	j=index_j(ig0)
979	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
980	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
981	ENDDO
982
983	IF (is_north_pole) then
984	DO i=1,iip1
985	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
986	ENDDO
987	ENDIF
988
989	IF (is_south_pole) then
990	DO i=1,iip1
991	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
992	ENDDO
993	ENDIF
994
995	ENDDO
996	c$OMP END DO NOWAIT
997	ENDDO
998
999	c 65. champ u:
1000	c ------------
1001	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1002	DO l=1,llm
1003	!CDIR ON_ADB(index_i)
1004	!CDIR ON_ADB(index_j)
1005	!cdir NODEP
1006	do ig0=kstart,kend
1007	i=index_i(ig0)
1008	j=index_j(ig0)
1009
1010	if (i/=iim) then
1011	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1012	endif
1013
1014	if (i==1) then
1015	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1016	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1017	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1018	endif
1019
1020	enddo
1021
1022	if (is_north_pole) then
1023	DO i=1,iip1
1024	pdufi(i,1,l) = 0.
1025	ENDDO
1026	endif
1027
1028	if (is_south_pole) then
1029	DO i=1,iip1
1030	pdufi(i,jjp1,l) = 0.
1031	ENDDO
1032	endif
1033
1034	ENDDO
1035	c$OMP END DO NOWAIT
1036
1037	c 67. champ v:
1038	c ------------
1039
1040	kstart=1
1041	kend=klon
1042
1043	if (is_north_pole) kstart=2
1044	if (is_south_pole) kend=klon-1-iim
1045
1046	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1047	DO l=1,llm
1048	!CDIR ON_ADB(index_i)
1049	!CDIR ON_ADB(index_j)
1050	!cdir NODEP
1051	do ig0=kstart,kend
1052	i=index_i(ig0)
1053	j=index_j(ig0)
1054	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1055	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1056	$ zdvfi2(ig0+iim,l))
1057	$ *cv(i,j)
1058	enddo
1059
1060	ENDDO
1061	c$OMP END DO NOWAIT
1062
1063
1064	c 68. champ v pres des poles:
1065	c ---------------------------
1066	c v = U * cos(long) + V * SIN(long)
1067
1068	if (is_north_pole) then
1069
1070	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1071	DO l=1,llm
1072
1073	DO i=1,iim
1074	pdvfi(i,1,l)=
1075	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1076
1077	pdvfi(i,1,l)=
1078	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1079	ENDDO
1080
1081	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1082
1083	ENDDO
1084	c$OMP END DO NOWAIT
1085
1086	endif
1087
1088	if (is_south_pole) then
1089
1090	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1091	DO l=1,llm
1092
1093	DO i=1,iim
1094	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1095	$ +zdvfi(klon,l)*SIN(rlonv(i))
1096
1097	pdvfi(i,jjm,l)=
1098	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1099	ENDDO
1100
1101	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1102
1103	ENDDO
1104	c$OMP END DO NOWAIT
1105
1106	endif
1107	c-----------------------------------------------------------------------
1108
1109	700 CONTINUE
1110
1111	firstcal = .FALSE.
1112
1113	#else
1114	write(lunout,*)
1115	& "calfis_p: for now can only work with parallel physics"
1116	stop
1117	#endif
1118	! of #ifdef CPP_PHYS
1119	#endif
1120	! of #ifdef CPP_PARA
1121	END

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

Download in other formats:

Original Format