Context Navigation

calfis_loc.F @ 1749

Last change on this file since 1749 was 1749, checked in by Ehouarn Millour, 12 years ago

Added handling of Newtonian (and Shallow Water) modes in dyn3dmem:

adapted calfis_loc.F and gcm.F
removed unused routines divgrad_p.F and gradiv_p.F
adapted iniacademic.F90 and sw_case_williamson.F to work in parallel and renamed them iniacademi_loc.F90 and sw_case_williamson_loc.F
fixed bug in exner_milieu_loc.F (filtreg_p form mod_filtreg_p should be used)

File size: 31.0 KB

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

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

Download in other formats:

Original Format