Context Navigation

calfis_loc.F @ 4464

Last change on this file since 4464 was 4464, checked in by lguez, 15 months ago
Replace stop by call to `abort_(physiq\|gcm)`
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:keywords set to `Id`
File size: 32.9 KB

Line
1	!
2	! $Id: calfis_loc.F 4464 2023-03-09 17:03:51Z lguez $
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	$ pdufi,
24	$ pdvfi,
25	$ pdhfi,
26	$ pdqfi,
27	$ pdpsfi)
28	#ifdef CPP_PHYS
29	! If using physics
30	c
31	c Auteur : P. Le Van, F. Hourdin
32	c .........
33	USE dimphy
34	USE mod_phys_lmdz_mpi_data, mpi_root_xx=>mpi_master
35	USE mod_phys_lmdz_omp_data, ONLY: klon_omp, klon_omp_begin
36	USE mod_const_mpi, ONLY: COMM_LMDZ
37	USE mod_interface_dyn_phys
38	USE IOPHY
39	#endif
40	#ifdef CPP_PARA
41	USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
42	$ ,jj_begin_dyn=>jj_begin,jj_end_dyn=>jj_end
43	USE Write_Field
44	Use Write_field_p
45	USE Times
46	#endif
47	USE infotrac, ONLY: nqtot, tracers
48	USE control_mod, ONLY: planet_type, nsplit_phys
49	#ifdef CPP_PHYS
50	USE callphysiq_mod, ONLY: call_physiq
51	#endif
52	USE comvert_mod, ONLY: preff, presnivs
53	USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
54
55	#ifdef CPP_PARA
56	IMPLICIT NONE
57	c=======================================================================
58	c
59	c 1. rearrangement des tableaux et transformation
60	c variables dynamiques > variables physiques
61	c 2. calcul des termes physiques
62	c 3. retransformation des tendances physiques en tendances dynamiques
63	c
64	c remarques:
65	c ----------
66	c
67	c - les vents sont donnes dans la physique par leurs composantes
68	c naturelles.
69	c - la variable thermodynamique de la physique est une variable
70	c intensive : T
71	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
72	c - les deux seules variables dependant de la geometrie necessaires
73	c pour la physique sont la latitude pour le rayonnement et
74	c l'aire de la maille quand on veut integrer une grandeur
75	c horizontalement.
76	c - les points de la physique sont les points scalaires de la
77	c la dynamique; numerotation:
78	c 1 pour le pole nord
79	c (jjm-1)*iim pour l'interieur du domaine
80	c ngridmx pour le pole sud
81	c ---> ngridmx=2+(jjm-1)*iim
82	c
83	c Input :
84	c -------
85	c ecritphy frequence d'ecriture (en jours)de histphy
86	c pucov covariant zonal velocity
87	c pvcov covariant meridional velocity
88	c pteta potential temperature
89	c pps surface pressure
90	c pmasse masse d'air dans chaque maille
91	c pts surface temperature (K)
92	c callrad clef d'appel au rayonnement
93	c
94	c Output :
95	c --------
96	c pdufi tendency for the natural zonal velocity (ms-1)
97	c pdvfi tendency for the natural meridional velocity
98	c pdhfi tendency for the potential temperature
99	c pdtsfi tendency for the surface temperature
100	c
101	c pdtrad radiative tendencies \ both input
102	c pfluxrad radiative fluxes / and output
103	c
104	c=======================================================================
105	c
106	c-----------------------------------------------------------------------
107	c
108	c 0. Declarations :
109	c ------------------
110
111	include "dimensions.h"
112	include "paramet.h"
113
114	INTEGER ngridmx
115	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
116
117	include "comgeom2.h"
118	include "iniprint.h"
119	#ifdef CPP_MPI
120	include 'mpif.h'
121	#endif
122	c Arguments :
123	c -----------
124	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
125	REAL,INTENT(IN):: jD_cur, jH_cur
126	REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
127	REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
128	REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
129	REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
130	REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
131	REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
132	REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential
133
134	REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
135	REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
136	REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
137	REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used
138
139	REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
140	REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
141	REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
142	REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
143
144	! tendencies (in */s) from the physics
145	REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
146	REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
147	REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
148	REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
149	REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)
150
151	#ifdef CPP_PHYS
152	! Ehouarn: for now calfis_p needs some informations from physics to compile
153	c Local variables :
154	c -----------------
155
156	INTEGER i,j,l,ig0,ig,iq,itr
157	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
158	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
159	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
160	c
161	REAL zrot(iip1,jjb_v:jje_v,llm) ! AdlC May 2014
162	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:), zrfi(:,:)
163	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
164	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
165	c
166	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
167	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
168	c
169	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
170	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
171	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
172	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
173
174	c
175	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
176	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
177	REAL,ALLOCATABLE,SAVE :: zpk_omp(:,:)
178	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
179	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
180	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
181	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
182	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
183	REAL,ALLOCATABLE,SAVE :: zrfi_omp(:,:)
184	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
185	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
186	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
187	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
188	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
189	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
190	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
191	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
192
193	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194	! Introduction du splitting (FH)
195	! Question pour Yann :
196	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
197	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
198	! soit allocatable (plutot par exemple que de passer une dimension
199	! dépendant du process en argument des routines) et que, du coup,
200	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
201	! Tu confirmes ?
202	! J'ai suivi le même principe pour les zdufic_omp
203	! Mais c'est surement bien que tu controles.
204	!
205
206	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
207	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
208	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
209	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
210	REAL jH_cur_split,zdt_split
211	LOGICAL debut_split,lafin_split
212	INTEGER isplit
213	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
214
215	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp,
216	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
217	c$OMP+ zrfi_omp,zqfi_omp,zdufi_omp,zdvfi_omp,
218	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
219	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
220
221	LOGICAL,SAVE :: first_omp=.true.
222	c$OMP THREADPRIVATE(first_omp)
223
224	REAL zsin(iim),zcos(iim),z1(iim)
225	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
226	REAL unskap, pksurcp
227	c
228	REAL SSUM
229
230	LOGICAL,SAVE :: firstcal=.true., debut=.true.
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	INTEGER :: jjb,jje
245
246	c
247	c-----------------------------------------------------------------------
248	c
249	c 1. Initialisations :
250	c --------------------
251	c
252
253	klon=klon_mpi
254
255	c
256	IF ( firstcal ) THEN
257	debut = .TRUE.
258	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
259	write(lunout,*) 'STOP dans calfis'
260	write(lunout,*)
261	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
262	write(lunout,*) ' ngridmx jjm iim '
263	write(lunout,*) ngridmx,jjm,iim
264	call abort_gcm("calfis_loc", "", 1)
265	ENDIF
266	c$OMP MASTER
267	ALLOCATE(zpsrf(klon))
268	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
269	ALLOCATE(zphi(klon,llm),zphis(klon))
270	ALLOCATE(zufi(klon,llm), zvfi(klon,llm),zrfi(klon,llm))
271	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
272	ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
273	ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
274	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
275	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
276	ALLOCATE(zdpsrf(klon))
277	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
278	ALLOCATE(flxwfi(klon,llm))
279	ALLOCATE(zpk(klon,llm))
280	c$OMP END MASTER
281	c$OMP BARRIER
282	ELSE
283	debut = .FALSE.
284	ENDIF
285
286	c
287	c
288	c-----------------------------------------------------------------------
289	c 40. transformation des variables dynamiques en variables physiques:
290	c ---------------------------------------------------------------
291
292	c 41. pressions au sol (en Pascals)
293	c ----------------------------------
294
295	c$OMP MASTER
296	call start_timer(timer_physic)
297	c$OMP END MASTER
298
299	c$OMP MASTER
300	!CDIR ON_ADB(index_i)
301	!CDIR ON_ADB(index_j)
302	do ig0=1,klon
303	i=index_i(ig0)
304	j=index_j(ig0)
305	zpsrf(ig0)=pps(i,j)
306	enddo
307	c$OMP END MASTER
308
309
310	c 42. pression intercouches :
311	c
312	c -----------------------------------------------------------------
313	c .... zplev definis aux (llm +1) interfaces des couches ....
314	c .... zplay definis aux ( llm ) milieux des couches ....
315	c -----------------------------------------------------------------
316
317	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
318	c
319	unskap = 1./ kappa
320	c
321	c print *,omp_rank,'klon--->',klon
322	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
323	DO l = 1, llmp1
324	!CDIR ON_ADB(index_i)
325	!CDIR ON_ADB(index_j)
326	do ig0=1,klon
327	i=index_i(ig0)
328	j=index_j(ig0)
329	zplev( ig0,l ) = pp(i,j,l)
330	enddo
331	ENDDO
332	c$OMP END DO NOWAIT
333
334	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
335	DO l=1,llm
336	do ig0=1,klon
337	i=index_i(ig0)
338	j=index_j(ig0)
339	zpk(ig0,l)=ppk(i,j,l)
340	enddo
341	ENDDO
342	c$OMP END DO NOWAIT
343
344	c
345	c
346
347	c 43. temperature naturelle (en K) et pressions milieux couches .
348	c ---------------------------------------------------------------
349	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
350	DO l=1,llm
351	!CDIR ON_ADB(index_i)
352	!CDIR ON_ADB(index_j)
353	do ig0=1,klon
354	i=index_i(ig0)
355	j=index_j(ig0)
356	pksurcp = ppk(i,j,l) / cpp
357	zplay(ig0,l) = preff * pksurcp ** unskap
358	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
359	enddo
360
361	ENDDO
362	c$OMP END DO NOWAIT
363
364	c 43.bis traceurs
365	c ---------------
366	c
367
368	itr = 0
369	DO iq=1,nqtot
370	IF(.NOT.tracers(iq)%isAdvected) CYCLE
371	itr = itr + 1
372	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
373	DO l=1,llm
374	!CDIR ON_ADB(index_i)
375	!CDIR ON_ADB(index_j)
376	do ig0=1,klon
377	i=index_i(ig0)
378	j=index_j(ig0)
379	zqfi(ig0,l,itr) = pq(i,j,l,iq)
380	enddo
381	ENDDO
382	c$OMP END DO NOWAIT
383	ENDDO
384
385
386	c Geopotentiel calcule par rapport a la surface locale:
387	c -----------------------------------------------------
388
389	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
390	DO l=1,llm
391	!CDIR ON_ADB(index_i)
392	!CDIR ON_ADB(index_j)
393	do ig0=1,klon
394	i=index_i(ig0)
395	j=index_j(ig0)
396	zphi(ig0,l) = pphi(i,j,l)
397	enddo
398	ENDDO
399	c$OMP END DO NOWAIT
400
401	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
402
403	c$OMP MASTER
404	!CDIR ON_ADB(index_i)
405	!CDIR ON_ADB(index_j)
406	do ig0=1,klon
407	i=index_i(ig0)
408	j=index_j(ig0)
409	zphis(ig0) = pphis(i,j)
410	enddo
411	c$OMP END MASTER
412
413
414	c CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
415
416	c$OMP BARRIER
417
418	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
419	DO l=1,llm
420	DO ig=1,klon
421	zphi(ig,l)=zphi(ig,l)-zphis(ig)
422	ENDDO
423	ENDDO
424	c$OMP END DO NOWAIT
425
426
427	c
428	c 45. champ u:
429	c ------------
430
431	kstart=1
432	kend=klon
433
434	if (is_north_pole_dyn) kstart=2
435	if (is_south_pole_dyn) kend=klon-1
436
437	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
438	DO l=1,llm
439	!CDIR ON_ADB(index_i)
440	!CDIR ON_ADB(index_j)
441	!CDIR SPARSE
442	do ig0=kstart,kend
443	i=index_i(ig0)
444	j=index_j(ig0)
445	if (i==1) then
446	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
447	$ + pucov(1,j,l)/cu(1,j) )
448	else
449	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
450	$ + pucov(i,j,l)/cu(i,j) )
451	endif
452	enddo
453	ENDDO
454	c$OMP END DO NOWAIT
455
456	c
457	C Alvaro de la Camara (May 2014)
458	C 46.1 Calcul de la vorticite et passage sur la grille physique
459	C --------------------------------------------------------------
460
461	jjb=jj_begin_dyn-1
462	jje=jj_end_dyn+1
463	if (is_north_pole_dyn) jjb=1
464	if (is_south_pole_dyn) jje=jjm
465
466	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
467
468	DO l=1,llm
469	do i=1,iim
470	do j=jjb,jje
471	zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
472	$ + pucov(i,j+1,l) - pucov(i,j,l))
473	$ / (cu(i,j)+cu(i,j+1))
474	$ / (cv(i+1,j)+cv(i,j)) *4
475	enddo
476	enddo
477	ENDDO
478
479
480	c 46.2champ v:
481	c -----------
482
483	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
484	DO l=1,llm
485	!CDIR ON_ADB(index_i)
486	!CDIR ON_ADB(index_j)
487	DO ig0=kstart,kend
488	i=index_i(ig0)
489	j=index_j(ig0)
490	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
491	$ + pvcov(i,j,l)/cv(i,j) )
492	if (j==1 .OR. j==jjp1) then ! AdlC MAY 2014
493	zrfi(ig0,l) = 0 ! AdlC MAY 2014
494	else
495	if(i==1)then
496	zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
497	$ +zrot(1,j-1,l)+zrot(1,j,l)) ! AdlC MAY 2014
498	else
499	zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
500	$ +zrot(i,j-1,l)+zrot(i,j,l)) ! AdlC MAY 2014
501	endif
502	endif
503
504
505	ENDDO
506	ENDDO
507	c$OMP END DO NOWAIT
508
509	c 47. champs de vents aux pole nord
510	c ------------------------------
511	c U = 1 / pi * integrale [ v * cos(long) * d long ]
512	c V = 1 / pi * integrale [ v * sin(long) * d long ]
513
514	if (is_north_pole_dyn) then
515	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
516	DO l=1,llm
517
518	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
519	DO i=2,iim
520	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
521	ENDDO
522
523	DO i=1,iim
524	zcos(i) = COS(rlonv(i))*z1(i)
525	zsin(i) = SIN(rlonv(i))*z1(i)
526	ENDDO
527
528	zufi(1,l) = SSUM(iim,zcos,1)/pi
529	zvfi(1,l) = SSUM(iim,zsin,1)/pi
530	zrfi(1,l) = 0.
531
532	ENDDO
533	c$OMP END DO NOWAIT
534	endif
535
536
537	c 48. champs de vents aux pole sud:
538	c ---------------------------------
539	c U = 1 / pi * integrale [ v * cos(long) * d long ]
540	c V = 1 / pi * integrale [ v * sin(long) * d long ]
541
542	if (is_south_pole_dyn) then
543	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
544	DO l=1,llm
545
546	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
547	DO i=2,iim
548	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
549	ENDDO
550
551	DO i=1,iim
552	zcos(i) = COS(rlonv(i))*z1(i)
553	zsin(i) = SIN(rlonv(i))*z1(i)
554	ENDDO
555
556	zufi(klon,l) = SSUM(iim,zcos,1)/pi
557	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
558	zrfi(klon,l) = 0.
559	ENDDO
560	c$OMP END DO NOWAIT
561	endif
562
563	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
564	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
565	DO l=1,llm
566	!CDIR ON_ADB(index_i)
567	!CDIR ON_ADB(index_j)
568	do ig0=1,klon
569	i=index_i(ig0)
570	j=index_j(ig0)
571	flxwfi(ig0,l) = flxw(i,j,l)
572	enddo
573	ENDDO
574	c$OMP END DO NOWAIT
575
576	c CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
577
578	c-----------------------------------------------------------------------
579	c Appel de la physique:
580	c ---------------------
581
582
583	c$OMP BARRIER
584	if (first_omp) then
585	klon=klon_omp
586
587	allocate(zplev_omp(klon,llm+1))
588	allocate(zplay_omp(klon,llm))
589	allocate(zpk_omp(klon,llm))
590	allocate(zphi_omp(klon,llm))
591	allocate(zphis_omp(klon))
592	allocate(presnivs_omp(llm))
593	allocate(zufi_omp(klon,llm))
594	allocate(zvfi_omp(klon,llm))
595	allocate(zrfi_omp(klon,llm)) ! LG Ari 2014
596	allocate(ztfi_omp(klon,llm))
597	allocate(zqfi_omp(klon,llm,nqtot))
598	allocate(zdufi_omp(klon,llm))
599	allocate(zdvfi_omp(klon,llm))
600	allocate(zdtfi_omp(klon,llm))
601	allocate(zdqfi_omp(klon,llm,nqtot))
602	allocate(zdufic_omp(klon,llm))
603	allocate(zdvfic_omp(klon,llm))
604	allocate(zdtfic_omp(klon,llm))
605	allocate(zdqfic_omp(klon,llm,nqtot))
606	allocate(zdpsrf_omp(klon))
607	allocate(flxwfi_omp(klon,llm))
608	first_omp=.false.
609	endif
610
611
612	klon=klon_omp
613	offset=klon_omp_begin-1
614
615	do l=1,llm+1
616	do i=1,klon
617	zplev_omp(i,l)=zplev(offset+i,l)
618	enddo
619	enddo
620
621	do l=1,llm
622	do i=1,klon
623	zplay_omp(i,l)=zplay(offset+i,l)
624	enddo
625	enddo
626
627	do l=1,llm
628	do i=1,klon
629	zpk_omp(i,l)=zpk(offset+i,l)
630	enddo
631	enddo
632
633	do l=1,llm
634	do i=1,klon
635	zphi_omp(i,l)=zphi(offset+i,l)
636	enddo
637	enddo
638
639	do i=1,klon
640	zphis_omp(i)=zphis(offset+i)
641	enddo
642
643
644	do l=1,llm
645	presnivs_omp(l)=presnivs(l)
646	enddo
647
648	do l=1,llm
649	do i=1,klon
650	zufi_omp(i,l)=zufi(offset+i,l)
651	enddo
652	enddo
653
654	do l=1,llm
655	do i=1,klon
656	zvfi_omp(i,l)=zvfi(offset+i,l)
657	enddo
658	enddo
659
660	do l=1,llm
661	do i=1,klon
662	zrfi_omp(i,l)=zrfi(offset+i,l)
663	enddo
664	enddo
665
666	do l=1,llm
667	do i=1,klon
668	ztfi_omp(i,l)=ztfi(offset+i,l)
669	enddo
670	enddo
671
672	do iq=1,nqtot
673	do l=1,llm
674	do i=1,klon
675	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
676	enddo
677	enddo
678	enddo
679
680	do l=1,llm
681	do i=1,klon
682	zdufi_omp(i,l)=zdufi(offset+i,l)
683	enddo
684	enddo
685
686	do l=1,llm
687	do i=1,klon
688	zdvfi_omp(i,l)=zdvfi(offset+i,l)
689	enddo
690	enddo
691
692	do l=1,llm
693	do i=1,klon
694	zdtfi_omp(i,l)=zdtfi(offset+i,l)
695	enddo
696	enddo
697
698	do iq=1,nqtot
699	do l=1,llm
700	do i=1,klon
701	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
702	enddo
703	enddo
704	enddo
705
706	do i=1,klon
707	zdpsrf_omp(i)=zdpsrf(offset+i)
708	enddo
709
710	do l=1,llm
711	do i=1,klon
712	flxwfi_omp(i,l)=flxwfi(offset+i,l)
713	enddo
714	enddo
715
716	c$OMP BARRIER
717
718
719	!$OMP MASTER
720	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
721	!$OMP END MASTER
722	zdt_split=dtphys/nsplit_phys
723	zdufic_omp(:,:)=0.
724	zdvfic_omp(:,:)=0.
725	zdtfic_omp(:,:)=0.
726	zdqfic_omp(:,:,:)=0.
727
728	#ifdef CPP_PHYS
729	do isplit=1,nsplit_phys
730
731	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
732	debut_split=debut.and.isplit==1
733	lafin_split=lafin.and.isplit==nsplit_phys
734
735	CALL call_physiq(klon,llm,nqtot,tracers(:)%name,
736	& debut_split,lafin_split,
737	& jD_cur,jH_cur_split,zdt_split,
738	& zplev_omp,zplay_omp,
739	& zpk_omp,zphi_omp,zphis_omp,
740	& presnivs_omp,
741	& zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
742	& flxwfi_omp,pducov,
743	& zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
744	& zdpsrf_omp)
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_dyn) 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_dyn) 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_dyn) kstart=2
969	if (is_south_pole_dyn) 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_dyn) then
982	DO i=1,iip1
983	pdpsfi(i,1) = zdpsrf(1)
984	enddo
985	endif
986
987	if (is_south_pole_dyn) 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_dyn) kstart=2
1003	if (is_south_pole_dyn) 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_dyn) 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_dyn) 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_dyn) 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_dyn) 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	itr = 0
1074	DO iq=1,nqtot
1075	IF(.NOT.tracers(iq)%isAdvected) CYCLE
1076	itr = itr + 1
1077	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1078	DO l=1,llm
1079	!CDIR ON_ADB(index_i)
1080	!CDIR ON_ADB(index_j)
1081	!cdir NODEP
1082	DO ig0=kstart,kend
1083	i=index_i(ig0)
1084	j=index_j(ig0)
1085	pdqfi(i,j,l,iq) = zdqfi(ig0,l,itr)
1086	if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,itr)
1087	ENDDO
1088
1089	IF (is_north_pole_dyn) then
1090	DO i=1,iip1
1091	pdqfi(i,1,l,iq) = zdqfi(1,l,itr)
1092	ENDDO
1093	ENDIF
1094
1095	IF (is_south_pole_dyn) then
1096	DO i=1,iip1
1097	pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,itr)
1098	ENDDO
1099	ENDIF
1100
1101	ENDDO
1102	c$OMP END DO NOWAIT
1103	ENDDO
1104
1105	c 65. champ u:
1106	c ------------
1107	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1108	DO l=1,llm
1109	!CDIR ON_ADB(index_i)
1110	!CDIR ON_ADB(index_j)
1111	!cdir NODEP
1112	do ig0=kstart,kend
1113	i=index_i(ig0)
1114	j=index_j(ig0)
1115
1116	if (i/=iim) then
1117	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1118	endif
1119
1120	if (i==1) then
1121	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
1122	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
1123	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
1124	endif
1125
1126	enddo
1127
1128	if (is_north_pole_dyn) then
1129	DO i=1,iip1
1130	pdufi(i,1,l) = 0.
1131	ENDDO
1132	endif
1133
1134	if (is_south_pole_dyn) then
1135	DO i=1,iip1
1136	pdufi(i,jjp1,l) = 0.
1137	ENDDO
1138	endif
1139
1140	ENDDO
1141	c$OMP END DO NOWAIT
1142
1143	c 67. champ v:
1144	c ------------
1145
1146	kstart=1
1147	kend=klon
1148
1149	if (is_north_pole_dyn) kstart=2
1150	if (is_south_pole_dyn) kend=klon-1-iim
1151
1152	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1153	DO l=1,llm
1154	!CDIR ON_ADB(index_i)
1155	!CDIR ON_ADB(index_j)
1156	!cdir NODEP
1157	do ig0=kstart,kend
1158	i=index_i(ig0)
1159	j=index_j(ig0)
1160	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
1161	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
1162	$ zdvfi2(ig0+iim,l))
1163	$ *cv(i,j)
1164	enddo
1165
1166	ENDDO
1167	c$OMP END DO NOWAIT
1168
1169
1170	c 68. champ v pres des poles:
1171	c ---------------------------
1172	c v = U * cos(long) + V * SIN(long)
1173
1174	if (is_north_pole_dyn) then
1175
1176	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1177	DO l=1,llm
1178
1179	DO i=1,iim
1180	pdvfi(i,1,l)=
1181	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
1182
1183	pdvfi(i,1,l)=
1184	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
1185	ENDDO
1186
1187	pdvfi(iip1,1,l) = pdvfi(1,1,l)
1188
1189	ENDDO
1190	c$OMP END DO NOWAIT
1191
1192	endif
1193
1194	if (is_south_pole_dyn) then
1195
1196	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1197	DO l=1,llm
1198
1199	DO i=1,iim
1200	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
1201	$ +zdvfi(klon,l)*SIN(rlonv(i))
1202
1203	pdvfi(i,jjm,l)=
1204	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
1205	ENDDO
1206
1207	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
1208
1209	ENDDO
1210	c$OMP END DO NOWAIT
1211
1212	endif
1213	c-----------------------------------------------------------------------
1214
1215	700 CONTINUE
1216
1217	firstcal = .FALSE.
1218
1219	#else
1220	call abort_gcm("calfis_loc",
1221	& "calfis_p: for now can only work with parallel physics", 1)
1222	#endif
1223	! of #ifdef CPP_PHYS
1224	#endif
1225	! of #ifdef CPP_PARA
1226	END

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

Download in other formats:

Original Format