Context Navigation

calfis_loc.F @ 1763

Last change on this file since 1763 was 1749, checked in by Ehouarn Millour, 12 years ago

Added handling of Newtonian (and Shallow Water) modes in dyn3dmem:

adapted calfis_loc.F and gcm.F
removed unused routines divgrad_p.F and gradiv_p.F
adapted iniacademic.F90 and sw_case_williamson.F to work in parallel and renamed them iniacademi_loc.F90 and sw_case_williamson_loc.F
fixed bug in exner_milieu_loc.F (filtreg_p form mod_filtreg_p should be used)

File size: 31.0 KB

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

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

Download in other formats:

Original Format