Context Navigation

calfis_loc.F @ 2272

Last change on this file since 2272 was 2258, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes 2216:2237 into testing branch
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: 31.1 KB

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

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

Download in other formats:

Original Format