Context Navigation

calfis_p.F @ 1972

Last change on this file since 1972 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.3 KB

Rev	Line
[630]	1	!
[1279]	2	! $Id: calfis_p.F 1907 2013-11-26 13:10:46Z fhourdin $
[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
[1146]	40	USE infotrac
[1403]	41	USE control_mod
[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 -----------
	114	LOGICAL lafin
[1403]	115	! REAL heure
	116	REAL, intent(in):: jD_cur, jH_cur
[630]	117	REAL pvcov(iip1,jjm,llm)
	118	REAL pucov(iip1,jjp1,llm)
	119	REAL pteta(iip1,jjp1,llm)
	120	REAL pmasse(iip1,jjp1,llm)
[1146]	121	REAL pq(iip1,jjp1,llm,nqtot)
[630]	122	REAL pphis(iip1,jjp1)
	123	REAL pphi(iip1,jjp1,llm)
	124	c
	125	REAL pdvcov(iip1,jjm,llm)
	126	REAL pducov(iip1,jjp1,llm)
	127	REAL pdteta(iip1,jjp1,llm)
[1146]	128	REAL pdq(iip1,jjp1,llm,nqtot)
[1403]	129	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
[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)
[1146]	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
[1617]	145	#ifdef CPP_PHYS
	146	! Ehouarn: for now calfis_p needs some informations from physics to compile
[630]	147	c Local variables :
	148	c -----------------
	149
	150	INTEGER i,j,l,ig0,ig,iq,iiq
[764]	151	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
	152	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
	153	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
[630]	154	c
[764]	155	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
	156	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
[630]	157	c
[764]	158	REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
	159	REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
[630]	160	c
[764]	161	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
	162	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
	163	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
[985]	164	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
	165
[630]	166	c
[764]	167	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
	168	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
	169	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
	170	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
	171	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
	172	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
	173	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
	174	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
	175	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
	176	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
	177	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
	178	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
	179	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
	180	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
[985]	181	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
[764]	182
[1403]	183	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	184	! Introduction du splitting (FH)
	185	! Question pour Yann :
	186	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
	187	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
	188	! soit allocatable (plutot par exemple que de passer une dimension
	189	! dépendant du process en argument des routines) et que, du coup,
	190	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
	191	! Tu confirmes ?
	192	! J'ai suivi le même principe pour les zdufic_omp
	193	! Mais c'est surement bien que tu controles.
	194	!
	195
	196	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
	197	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
	198	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
	199	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
	200	REAL jH_cur_split,zdt_split
	201	LOGICAL debut_split,lafin_split
	202	INTEGER isplit
	203	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	204
[764]	205	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
	206	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
[985]	207	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
[1403]	208	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
	209	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
[764]	210
	211	LOGICAL,SAVE :: first_omp=.true.
	212	c$OMP THREADPRIVATE(first_omp)
	213
[630]	214	REAL zsin(iim),zcos(iim),z1(iim)
	215	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	216	REAL unskap, pksurcp
[764]	217	c
	218	cIM diagnostique PVteta, Amip2
	219	INTEGER ntetaSTD
	220	PARAMETER(ntetaSTD=3)
	221	REAL rtetaSTD(ntetaSTD)
[1617]	222	DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
[764]	223	REAL PVteta(klon,ntetaSTD)
	224
[630]	225
	226	REAL SSUM
	227
	228	LOGICAL firstcal, debut
	229	DATA firstcal/.true./
	230	SAVE firstcal,debut
[764]	231	c$OMP THREADPRIVATE(firstcal,debut)
[630]	232
[764]	233	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
[630]	234	INTEGER :: ierr
[1000]	235	#ifdef CPP_MPI
[630]	236	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
[1000]	237	#else
	238	INTEGER,dimension(1,4) :: Status
	239	#endif
[630]	240	INTEGER, dimension(4) :: Req
[764]	241	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
	242	integer :: k,kstart,kend
	243	INTEGER :: offset
[1654]	244
	245	LOGICAL tracerdyn
[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
[1615]	491	IF (is_sequential.and.(planet_type=="earth")) THEN
[1617]	492	#ifdef CPP_PHYS
	493	! PVtheta calls tetalevel, which is in the physics
[764]	494	cIM calcul PV a teta=350, 380, 405K
	495	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
	496	$ ztfi,zplay,zplev,
	497	$ ntetaSTD,rtetaSTD,PVteta)
	498	c
[1617]	499	#endif
[764]	500	ENDIF
[960]	501
	502	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
[630]	503	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
	504
	505	c-----------------------------------------------------------------------
	506	c Appel de la physique:
	507	c ---------------------
	508
	509
[764]	510	c$OMP BARRIER
	511	if (first_omp) then
[774]	512	klon=klon_omp
[764]	513
	514	allocate(zplev_omp(klon,llm+1))
	515	allocate(zplay_omp(klon,llm))
	516	allocate(zphi_omp(klon,llm))
	517	allocate(zphis_omp(klon))
	518	allocate(presnivs_omp(llm))
	519	allocate(zufi_omp(klon,llm))
	520	allocate(zvfi_omp(klon,llm))
	521	allocate(ztfi_omp(klon,llm))
[1146]	522	allocate(zqfi_omp(klon,llm,nqtot))
[764]	523	allocate(zdufi_omp(klon,llm))
	524	allocate(zdvfi_omp(klon,llm))
	525	allocate(zdtfi_omp(klon,llm))
[1146]	526	allocate(zdqfi_omp(klon,llm,nqtot))
[1403]	527	allocate(zdufic_omp(klon,llm))
	528	allocate(zdvfic_omp(klon,llm))
	529	allocate(zdtfic_omp(klon,llm))
	530	allocate(zdqfic_omp(klon,llm,nqtot))
[764]	531	allocate(zdpsrf_omp(klon))
[985]	532	allocate(flxwfi_omp(klon,llm))
[764]	533	first_omp=.false.
	534	endif
	535
	536
[774]	537	klon=klon_omp
	538	offset=klon_omp_begin-1
[764]	539
	540	do l=1,llm+1
	541	do i=1,klon
	542	zplev_omp(i,l)=zplev(offset+i,l)
	543	enddo
	544	enddo
	545
	546	do l=1,llm
	547	do i=1,klon
	548	zplay_omp(i,l)=zplay(offset+i,l)
	549	enddo
	550	enddo
	551
	552	do l=1,llm
	553	do i=1,klon
	554	zphi_omp(i,l)=zphi(offset+i,l)
	555	enddo
	556	enddo
	557
	558	do i=1,klon
	559	zphis_omp(i)=zphis(offset+i)
	560	enddo
	561
	562
	563	do l=1,llm
	564	presnivs_omp(l)=presnivs(l)
	565	enddo
	566
	567	do l=1,llm
	568	do i=1,klon
	569	zufi_omp(i,l)=zufi(offset+i,l)
	570	enddo
	571	enddo
	572
	573	do l=1,llm
	574	do i=1,klon
	575	zvfi_omp(i,l)=zvfi(offset+i,l)
	576	enddo
	577	enddo
	578
	579	do l=1,llm
	580	do i=1,klon
	581	ztfi_omp(i,l)=ztfi(offset+i,l)
	582	enddo
	583	enddo
	584
[1146]	585	do iq=1,nqtot
[764]	586	do l=1,llm
	587	do i=1,klon
	588	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
	589	enddo
	590	enddo
	591	enddo
	592
	593	do l=1,llm
	594	do i=1,klon
	595	zdufi_omp(i,l)=zdufi(offset+i,l)
	596	enddo
	597	enddo
	598
	599	do l=1,llm
	600	do i=1,klon
	601	zdvfi_omp(i,l)=zdvfi(offset+i,l)
	602	enddo
	603	enddo
	604
	605	do l=1,llm
	606	do i=1,klon
	607	zdtfi_omp(i,l)=zdtfi(offset+i,l)
	608	enddo
	609	enddo
	610
[1146]	611	do iq=1,nqtot
[764]	612	do l=1,llm
	613	do i=1,klon
	614	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
	615	enddo
	616	enddo
	617	enddo
	618
	619	do i=1,klon
	620	zdpsrf_omp(i)=zdpsrf(offset+i)
	621	enddo
[985]	622
	623	do l=1,llm
	624	do i=1,klon
	625	flxwfi_omp(i,l)=flxwfi(offset+i,l)
	626	enddo
	627	enddo
[764]	628
	629	c$OMP BARRIER
	630
[1403]	631	!$OMP MASTER
[1407]	632	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
[1403]	633	!$OMP END MASTER
	634	zdt_split=dtphys/nsplit_phys
	635	zdufic_omp(:,:)=0.
	636	zdvfic_omp(:,:)=0.
	637	zdtfic_omp(:,:)=0.
	638	zdqfic_omp(:,:,:)=0.
	639
[1615]	640	#ifdef CPP_PHYS
[1403]	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
[1654]	647	if (planet_type=="earth") then
[1403]	648
[630]	649	CALL physiq (klon,
	650	. llm,
[1403]	651	. debut_split,
	652	. lafin_split,
[1279]	653	. jD_cur,
[1403]	654	. jH_cur_split,
	655	. zdt_split,
[764]	656	. zplev_omp,
	657	. zplay_omp,
	658	. zphi_omp,
	659	. zphis_omp,
	660	. presnivs_omp,
[630]	661	. clesphy0,
[764]	662	. zufi_omp,
	663	. zvfi_omp,
	664	. ztfi_omp,
	665	. zqfi_omp,
[960]	666	c#ifdef INCA
[985]	667	. flxwfi_omp,
[960]	668	c#endif
[764]	669	. zdufi_omp,
	670	. zdvfi_omp,
	671	. zdtfi_omp,
	672	. zdqfi_omp,
	673	. zdpsrf_omp,
	674	cIM diagnostique PVteta, Amip2
	675	. pducov,
	676	. PVteta)
[1403]	677
[1654]	678	else if ( planet_type=="generic" ) then
	679
	680	CALL physiq (klon, !! ngrid
	681	. llm, !! nlayer
	682	. nqtot, !! nq
	683	. tname, !! tracer names from dynamical core (given in infotrac)
	684	. debut_split, !! firstcall
	685	. lafin_split, !! lastcall
[1676]	686	. jD_cur, !! pday. see leapfrog_p
[1654]	687	. jH_cur_split, !! ptime "fraction of day"
	688	. zdt_split, !! ptimestep
	689	. zplev_omp, !! pplev
	690	. zplay_omp, !! pplay
	691	. zphi_omp, !! pphi
	692	. zufi_omp, !! pu
	693	. zvfi_omp, !! pv
	694	. ztfi_omp, !! pt
	695	. zqfi_omp, !! pq
	696	. flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
	697	. zdufi_omp, !! pdu
	698	. zdvfi_omp, !! pdv
	699	. zdtfi_omp, !! pdt
	700	. zdqfi_omp, !! pdq
	701	. zdpsrf_omp, !! pdpsrf
	702	. tracerdyn) !! tracerdyn <-- utilite ???
	703
	704	endif ! of if (planet_type=="earth")
	705
	706
[1403]	707	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
	708	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
	709	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
	710	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
	711
	712	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
	713	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
	714	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
	715	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
	716
	717	enddo
	718
[1615]	719	#endif
	720	! of #ifdef CPP_PHYS
	721
[1403]	722	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
	723	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
	724	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
	725	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
[764]	726	c$OMP BARRIER
	727
	728	do l=1,llm+1
	729	do i=1,klon
	730	zplev(offset+i,l)=zplev_omp(i,l)
	731	enddo
	732	enddo
	733
	734	do l=1,llm
	735	do i=1,klon
	736	zplay(offset+i,l)=zplay_omp(i,l)
	737	enddo
	738	enddo
	739
	740	do l=1,llm
	741	do i=1,klon
	742	zphi(offset+i,l)=zphi_omp(i,l)
	743	enddo
	744	enddo
	745
	746
	747	do i=1,klon
	748	zphis(offset+i)=zphis_omp(i)
	749	enddo
	750
	751
	752	do l=1,llm
	753	presnivs(l)=presnivs_omp(l)
	754	enddo
	755
	756	do l=1,llm
	757	do i=1,klon
	758	zufi(offset+i,l)=zufi_omp(i,l)
	759	enddo
	760	enddo
	761
	762	do l=1,llm
	763	do i=1,klon
	764	zvfi(offset+i,l)=zvfi_omp(i,l)
	765	enddo
	766	enddo
	767
	768	do l=1,llm
	769	do i=1,klon
	770	ztfi(offset+i,l)=ztfi_omp(i,l)
	771	enddo
	772	enddo
	773
[1146]	774	do iq=1,nqtot
[764]	775	do l=1,llm
	776	do i=1,klon
	777	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
	778	enddo
	779	enddo
	780	enddo
	781
	782	do l=1,llm
	783	do i=1,klon
	784	zdufi(offset+i,l)=zdufi_omp(i,l)
	785	enddo
	786	enddo
	787
	788	do l=1,llm
	789	do i=1,klon
	790	zdvfi(offset+i,l)=zdvfi_omp(i,l)
	791	enddo
	792	enddo
	793
	794	do l=1,llm
	795	do i=1,klon
	796	zdtfi(offset+i,l)=zdtfi_omp(i,l)
	797	enddo
	798	enddo
	799
[1146]	800	do iq=1,nqtot
[764]	801	do l=1,llm
	802	do i=1,klon
	803	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
	804	enddo
	805	enddo
	806	enddo
	807
	808	do i=1,klon
	809	zdpsrf(offset+i)=zdpsrf_omp(i)
	810	enddo
	811
	812
	813	klon=klon_mpi
[630]	814	500 CONTINUE
[764]	815	c$OMP BARRIER
[630]	816
[764]	817	c$OMP MASTER
[630]	818	call stop_timer(timer_physic)
[764]	819	c$OMP END MASTER
[1000]	820
	821	IF (using_mpi) THEN
	822
[630]	823	if (MPI_rank>0) then
[764]	824
	825	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	826	DO l=1,llm
	827	du_send(1:iim,l)=zdufi(1:iim,l)
	828	dv_send(1:iim,l)=zdvfi(1:iim,l)
	829	ENDDO
	830	c$OMP END DO NOWAIT
	831
	832	c$OMP BARRIER
[1000]	833	#ifdef CPP_MPI
[764]	834	c$OMP MASTER
[985]	835	!$OMP CRITICAL (MPI)
[630]	836	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
[764]	837	& COMM_LMDZ,Req(1),ierr)
[630]	838	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
[764]	839	& COMM_LMDZ,Req(2),ierr)
[985]	840	!$OMP END CRITICAL (MPI)
[764]	841	c$OMP END MASTER
[1000]	842	#endif
[764]	843	c$OMP BARRIER
[630]	844
	845	endif
	846
	847	if (MPI_rank<MPI_Size-1) then
[764]	848	c$OMP BARRIER
[1000]	849	#ifdef CPP_MPI
[764]	850	c$OMP MASTER
[985]	851	!$OMP CRITICAL (MPI)
[630]	852	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
[764]	853	& COMM_LMDZ,Req(3),ierr)
[630]	854	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
[764]	855	& COMM_LMDZ,Req(4),ierr)
[985]	856	!$OMP END CRITICAL (MPI)
[764]	857	c$OMP END MASTER
[1000]	858	#endif
[630]	859	endif
[764]	860
	861	c$OMP BARRIER
[1000]	862
	863
	864	#ifdef CPP_MPI
[764]	865	c$OMP MASTER
[985]	866	!$OMP CRITICAL (MPI)
[630]	867	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
	868	call MPI_WAITALL(4,Req(1),Status,ierr)
	869	else if (MPI_rank>0) then
	870	call MPI_WAITALL(2,Req(1),Status,ierr)
	871	else if (MPI_rank <MPI_Size-1) then
	872	call MPI_WAITALL(2,Req(3),Status,ierr)
	873	endif
[985]	874	!$OMP END CRITICAL (MPI)
[764]	875	c$OMP END MASTER
[1000]	876	#endif
	877
[764]	878	c$OMP BARRIER
	879
[1000]	880	ENDIF ! using_mpi
	881
	882
[764]	883	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	884	DO l=1,llm
	885
	886	zdufi2(1:klon,l)=zdufi(1:klon,l)
	887	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
	888
	889	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
	890	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
	891
[774]	892	pdhfi(:,jj_begin,l)=0
	893	pdqfi(:,jj_begin,l,:)=0
	894	pdufi(:,jj_begin,l)=0
	895	pdvfi(:,jj_begin,l)=0
[764]	896
[774]	897	if (.not. is_south_pole) then
	898	pdhfi(:,jj_end,l)=0
	899	pdqfi(:,jj_end,l,:)=0
	900	pdufi(:,jj_end,l)=0
	901	pdvfi(:,jj_end,l)=0
[764]	902	endif
[630]	903
[764]	904	ENDDO
	905	c$OMP END DO NOWAIT
[630]	906
[764]	907	c$OMP MASTER
[774]	908	pdpsfi(:,jj_begin)=0
	909	if (.not. is_south_pole) then
	910	pdpsfi(:,jj_end)=0
[630]	911	endif
[764]	912	c$OMP END MASTER
[630]	913	c-----------------------------------------------------------------------
	914	c transformation des tendances physiques en tendances dynamiques:
	915	c ---------------------------------------------------------------
	916
	917	c tendance sur la pression :
	918	c -----------------------------------
	919	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
	920	c
	921	c 62. enthalpie potentielle
	922	c ---------------------
	923
	924	kstart=1
	925	kend=klon
	926
[774]	927	if (is_north_pole) kstart=2
	928	if (is_south_pole) kend=klon-1
[630]	929
[764]	930	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	931	DO l=1,llm
	932
[1279]	933	!CDIR ON_ADB(index_i)
	934	!CDIR ON_ADB(index_j)
	935	!cdir NODEP
[630]	936	do ig0=kstart,kend
[774]	937	i=index_i(ig0)
	938	j=index_j(ig0)
[630]	939	pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	940	if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	941	enddo
	942
[774]	943	if (is_north_pole) then
[630]	944	DO i=1,iip1
	945	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	946	enddo
	947	endif
	948
[774]	949	if (is_south_pole) then
[630]	950	DO i=1,iip1
	951	pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
	952	ENDDO
	953	endif
	954	ENDDO
[764]	955	c$OMP END DO NOWAIT
[630]	956
	957	c 62. humidite specifique
	958	c ---------------------
[1279]	959	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	960	! DO iq=1,nqtot
	961	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	962	! DO l=1,llm
	963	!!!cdir NODEP
	964	! do ig0=kstart,kend
	965	! i=index_i(ig0)
	966	! j=index_j(ig0)
	967	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
	968	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
	969	! enddo
	970	!
	971	! if (is_north_pole) then
	972	! do i=1,iip1
	973	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	974	! enddo
	975	! endif
	976	!
	977	! if (is_south_pole) then
	978	! do i=1,iip1
	979	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
	980	! enddo
	981	! endif
	982	! ENDDO
	983	!c$OMP END DO NOWAIT
	984	! ENDDO
[630]	985
	986	c 63. traceurs
	987	c ------------
	988	C initialisation des tendances
[764]	989
	990	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	991	DO l=1,llm
	992	pdqfi(:,:,l,:)=0.
	993	ENDDO
	994	c$OMP END DO NOWAIT
	995
[630]	996	C
[1279]	997	!cdir NODEP
[1146]	998	DO iq=1,nqtot
[630]	999	iiq=niadv(iq)
[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	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
	1009	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
	1010	ENDDO
	1011
[774]	1012	IF (is_north_pole) then
[630]	1013	DO i=1,iip1
	1014	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
	1015	ENDDO
	1016	ENDIF
	1017
[774]	1018	IF (is_south_pole) then
[630]	1019	DO i=1,iip1
	1020	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
	1021	ENDDO
	1022	ENDIF
	1023
	1024	ENDDO
[764]	1025	c$OMP END DO NOWAIT
[630]	1026	ENDDO
	1027
	1028	c 65. champ u:
	1029	c ------------
[764]	1030	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1031	DO l=1,llm
[1279]	1032	!CDIR ON_ADB(index_i)
	1033	!CDIR ON_ADB(index_j)
	1034	!cdir NODEP
[630]	1035	do ig0=kstart,kend
[774]	1036	i=index_i(ig0)
	1037	j=index_j(ig0)
[630]	1038
	1039	if (i/=iim) then
	1040	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
	1041	endif
	1042
	1043	if (i==1) then
	1044	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
	1045	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
[764]	1046	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
[630]	1047	endif
	1048
	1049	enddo
	1050
[774]	1051	if (is_north_pole) then
[630]	1052	DO i=1,iip1
	1053	pdufi(i,1,l) = 0.
	1054	ENDDO
	1055	endif
	1056
[774]	1057	if (is_south_pole) then
[630]	1058	DO i=1,iip1
	1059	pdufi(i,jjp1,l) = 0.
	1060	ENDDO
	1061	endif
	1062
	1063	ENDDO
[764]	1064	c$OMP END DO NOWAIT
[630]	1065
	1066	c 67. champ v:
	1067	c ------------
	1068
	1069	kstart=1
	1070	kend=klon
	1071
[774]	1072	if (is_north_pole) kstart=2
	1073	if (is_south_pole) kend=klon-1-iim
[630]	1074
[764]	1075	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1076	DO l=1,llm
[1279]	1077	!CDIR ON_ADB(index_i)
	1078	!CDIR ON_ADB(index_j)
	1079	!cdir NODEP
[630]	1080	do ig0=kstart,kend
[774]	1081	i=index_i(ig0)
	1082	j=index_j(ig0)
[630]	1083	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
	1084	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
	1085	$ zdvfi2(ig0+iim,l))
	1086	$ *cv(i,j)
	1087	enddo
	1088
	1089	ENDDO
[764]	1090	c$OMP END DO NOWAIT
[630]	1091
	1092
	1093	c 68. champ v pres des poles:
	1094	c ---------------------------
	1095	c v = U * cos(long) + V * SIN(long)
	1096
[774]	1097	if (is_north_pole) then
[764]	1098
	1099	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1100	DO l=1,llm
	1101
	1102	DO i=1,iim
	1103	pdvfi(i,1,l)=
	1104	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	1105
	1106	pdvfi(i,1,l)=
	1107	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	1108	ENDDO
	1109
	1110	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	1111
	1112	ENDDO
[764]	1113	c$OMP END DO NOWAIT
[630]	1114
	1115	endif
	1116
[774]	1117	if (is_south_pole) then
[764]	1118
	1119	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1120	DO l=1,llm
[630]	1121
	1122	DO i=1,iim
	1123	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
	1124	$ +zdvfi(klon,l)*SIN(rlonv(i))
	1125
	1126	pdvfi(i,jjm,l)=
	1127	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
	1128	ENDDO
	1129
	1130	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	1131
	1132	ENDDO
[764]	1133	c$OMP END DO NOWAIT
[630]	1134
	1135	endif
	1136	c-----------------------------------------------------------------------
	1137
	1138	700 CONTINUE
	1139
	1140	firstcal = .FALSE.
	1141
[1617]	1142	#else
	1143	write(lunout,*)
	1144	& "calfis_p: for now can only work with parallel physics"
	1145	stop
	1146	#endif
	1147	! of #ifdef CPP_PHYS
[630]	1148	RETURN
	1149	END

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

Download in other formats:

Original Format