Context Navigation

calfis_p.F @ 1249

Last change on this file since 1249 was 1231, checked in by lguez, 15 years ago

These are changes related to ozone in radiative transfer computations:

-- Moved "max(ozonecm, 1e-12)" from procedure "radlwsw" to procedure

"ozonecm".

-- Removed ratio "p_surface / p_reference" in the computation of

"pozon" in "radlwsw".

-- Removed ratio "p_reference / p_surface" in the computation of

"zoz" in "SW_LMDAR4" and "SW_AEROAR4".

-- Removed ratio "mass of air / mass of ozone" in the computation of

"zoz" in "LW_LMDAR4".

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 28.4 KB

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

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

Download in other formats:

Original Format