Context Navigation

calfis_p.F @ 873

Last change on this file since 873 was 849, checked in by emillour, 13 years ago

Work on common dynamics and interfacing with different physics:

Put calls to PVtheta in dynamics between CPP_EARTH flags (because it calls tetalevel, which is supposed to be in the physics; only OK for Earth...).
Adapted makelmdz script so that one can compile main programs in dyn* or phy* (makelmdz_fcm already capable of doing that).
Moved start2archive-VT.F and start2archive-VT.F to phyvenus (as start2archive.F and newstart.F); leave adapting them to Titan for later.
Small correction to phyvenus/testphys1d.F (use module control_mod instead of control.h and call disvert_noterre).
removed "use histcom" in phyvenus/physiq.F and phytitan/physiq.F ; it is not needed since there is already a "use ioipsl" (and it moreover confused makelmdz_fcm...)
Had to add declaration of variable "zlsm1" in phytitan/physiq.F because it is used in "write_hist.h"; but note that it is used while not initialized (but what should it be initialized to?).

File size: 31.4 KB

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

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

Download in other formats:

Original Format