Context Navigation

calfis_p.F @ 1329

Last change on this file since 1329 was 1325, checked in by idelkadi, 15 years ago
Decoupage du pas de temps physique en N sous-pas (splitting)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.8 KB

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

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

Download in other formats:

Original Format