Context Navigation

calfis_p.F @ 1720

Last change on this file since 1720 was 1407, checked in by Laurent Fairhead, 14 years ago

Some extraneous prints deleted for optimisation purposes

Des prints supprimés pour optimisation

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

File size: 28.0 KB

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

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

Download in other formats:

Original Format