Context Navigation

calfis_p.F @ 2037

Last change on this file since 2037 was 2037, checked in by lguez, 10 years ago

PVteta computed by PVtheta was not used. Also there were a couple of
problems with PVtheta:

-- PVtheta calls tetalevel in phylmd, and interpolates at

Earth-specific values of potential temperature.

-- PVtheta calls tourabs, and the computation of rot in tourabs should

be modified (it is not correct when there is a zoom).

-- Even when there is no zoom, the computation of rot in tourabs

should probably be modified: directly combine vcov and ucov, as in
tourpot, instead of dividing by cv and cu.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 29.7 KB

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

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

Download in other formats:

Original Format