Context Navigation

calfis_loc.F @ 2343

Last change on this file since 2343 was 2333, checked in by lguez, 9 years ago

New parameterization of gravity wave drag due to front/jet systems, by

de la Camara and F. Lott. The new Camara-Lott parameterization

replaces the Hines parameterization so it is activated if not ok_hines
and ok_gwd_rando.

Also changed distribution of phase speeds in FLOTT_GWD_rando, from
uniform to Gaussian. Bug fix in sugwd_strato. Bug fix in the arguments
of the call to add_phys_tend for methane oxydation.

For the new Camara-Lott parameterization, we need to compute relative
vorticity in calfis and pass it as a new argument "rot" to
physiq. Interpolation of relative vorticity to the physics grid is not
optimal for now: it is not weighted by cell areas.

Alvaro de la Camara, Fran\c{}cois Lott

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

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

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

Download in other formats:

Original Format