Context Navigation

calfis_p.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.3 KB

Line
1	!
2	! $Id: calfis_p.F 1907 2013-11-26 13:10:46Z lguez $
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
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